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ABSTRACT 


A  study  has  been  made  of  the  problem  of  estimating 
parameters  or  coefficients  in  systems  of  ordinary  differ¬ 
ential  equations  from  experimental  data.  The  study  was 
limited  to  a  special  class  of  problems.  However,  the  basic 
concepts  which  were  developed  can,  in  principle,  be  extended 
to  more  general  problems  without  much  difficulty. 

It  was  assumed  that  for  a  given  physical  system, 
several  experiments  have  been  performed  with  different,  but 
known  initial  conditions.  It  was  further  assumed  that  some 
or  all  of  the  state  variables  had  been  measured  at  one  or 
more  values  of  the  independent  variable.  Of  course,  it  is 
also  necessary  to  assume  that  the  given  physical  system  can 
be  adequately  modeled  by  a  system  of  ordinary  differential 
equations  which  is  known  to  within  a  few  parameters. 

This  identification  problem  was  formulated  as  a 
multipoint  boundary  value  problem  which  was  inevitably  non¬ 
linear.  Attempts  to  use  quasilinearization  directly  did  not 
always  converge.  To  alleviate  this  difficulty,  the  concept 
of  boundary  value  or  data  pertubation  was  introduced.  This 
concept  is  based  on  the  assumption  that  quasilinearization 
will  converge  provided  the  initial  guess  is  sufficiently 
close.  If  this  is  so,  then  pseudo-boundary  conditions  can  be 
set  up  near  the  trajectories  produced  by  the  initial  guess. 
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These  pseudo-boundary  conditions  are  chosen  so  that  the 
actual  solution  may  be  approached  in  a  finite  number  of 
steps,  each  step  being  solved  by  quasilinearization.  Thi 
procedure  was  found  satisfactory  for  a  number  of  examples 
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I.  INTRODUCTION 


In  recent  years,  with  the  exponential  advances  in 
technology,  it  has  been  increasingly  evident  that  more 
sophisticated  mathematical  models  must  be  employed  if 
industrial  processes  are  to  be  designed  and  operated  more 
efficiently. 

Management,  in  the  not  too  distant  past,  took  pride 
in  the  amount  of  "over  production"  that  was  obtained  from 
a  particular  plant.  However^  now  with  their  thinking  on  a 
more  sophisticated  level,  they  realize  that  the  "over 
production"  is  a  direct  result  of  inadequate  design  proce¬ 
dures  which  resulted  in  a  considerable  amount  of  investment 
capital  being  wasted  on  large  safety  factors. 

Recent  tendencies  toward  automation  are  completely 
dependent  on  having  adequate  mathematical  models.  A  large 
part  of  optimization  theory  assumes  that  the  models  exist; 
hence  the  development  of  techniques  for  obtaining  these 
models  is  of  utmost  importance. 

Finding  mathematical  models  which  are  adequate  for 
these  various  uses  is  not  a  trivial  problem.  A  purely 
theoretical  approach  is  sufficient  in  only  the  simplest  of 
cases.  In  more  complex  cases,  either  the  system  is  not  well 
enough  understood  or  its  model  is  so  complicated  that  it  can 
not  be  solved  in  a  reasonable  amount  of  time,  if  at  all.  The 
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definition  of  "a  reasonable  amount  of  time",  of  course  depends 
on  the  particular  problem.  With  these  complex  systems,  a 
correlative  procedure  is  needed. 

For  most  control  systems,  the  simplest  acceptable 
correlative  model  is  a  system  of  ordinary  differential  equa¬ 
tions.  Theoretical  considerations  usually  can  be  used  to 
derive  the  system  of  ordinary  differential  eauations  in  which 
only  a  few  parameters  need  to  be  determined. 

Similarly,  many  design  models  are  also  in  the  form  of 
a  system  of  ordinary  differential  equations,  known  except  for 
a  few  constant  parameters.  Plug  flow  and  batch  reactors  are 
the  examples  of  this  type  of  system;  these  are  used  extensive¬ 
ly  in  this  work. 

In  general,  data  taken  from  the  actual  system  or  a 
laboratory  scaled  model  of  the  system  are  used  to  identify  the 
unknown  parameters  in  the  mathematical  model . 

Unfortunately,  once  the  data  are  taken,  there  are  no 
completely  satisfactory  general  techniques  for  finding  these 
parameters.  The  techniques  that  do  exist  often  depend  on  some 
specific  characteristic  of  the  problem  (such  as  linearity) . 

A  good  summary  of  the  existing  techniques  with  an 
extensive  list  of  pertinent  references  has  been  supplied  by 
Cuenod  and  Gage  (5)  . 

Many  of  the  existing  techniques  require  that  the  experi¬ 
ments  be  performed  and  the  data  collected  in  a  predetermined 
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manner.  More  specifically,  the  procedure  for  identifying 
certain  control  systems  is  based  on  observed  deviations  of 
the  system  from  a  target  value  or  set  point  caused  by 
disturbances  introduced  into  the  system.  The  nature  of 
these  disturbances  is  dependent  on  the  particular  method  being 
employed.  An  example  is  the  sinusoidal  disturbance  required 
for  the  frequency  response  method. 

With  many  systems,  it  is  not  convenient  to  introduce 
arbitrary  disturbances  which  are  a  function  of  the  independent 
variable.  This  is  especially  true  with  the  identification  of 
certain  design  models. 

There  are  a  few  methods  for  handling  the  identification 
problem  which  are  not  as  dependent  on  the  types  of  disturbances 
employed,  provided  that  they  adequately  excite  all  the  different 
modes  of  the  system.  To  discuss  some  of  these  methods,  a  system 
will  be  considered  described  below 


dy 

—  =  g  (y,  a,  t)  (1-1) 

dt 


with  y  being  the  state  vector,  t  the  independent  variable,  and 

a  the  unknown  parameter  vector . 

One  approach  that  is  quite  favored  because  of  its 

dy 

simplicity  is  to  measure  y  and  —  as  functions  of  the 

dt 

independent  variable  and  minimize  the  objective  function 


tf 


o 


dy 

dt 


-  g;  (y,  a,  t) 


z 


dt 


(1-2) 
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li  M  •  dZ 

with  1 1 x| |  being  the  length  of  the  vector  x,  —  being  the  measured 

*  dt 

value  of  the  derivative  and  y_  being  the  measured  value  of  the 
state  vector.  If  the  measurements  are  not  taken  continuously, 
but  at  discrete  values  of  the  independent  variable,  the  cor¬ 
responding  objective  function  is 

A 

dy_.  .  2 

z  =  ?  U  —  -  2.  <Zi>  2  V  I  I  (1-3) 

i  dt 


The  principal  objection  to  this  method  is  that  it  is  quite 
sensitive  to  noise  in  the  data  which  make  it  difficult  to 

A 

d  y_ 

determine  accurate  values  of  —  .  For  some  systems,  the 

dy 

direct  measurement  of  —  is  not  difficult  and  this  method 

dt 

can  be  employed  profitably. 

An  alternative  approach  is  to  reformulate  the  problem 
by  letting  the  vector  a  be  a  function  of  the  independent 
variable  which  is  a  solution  to  the  differential  equation 


da 

—  =0  (1-4) 

dt 

Doing  this,  the  identification  problem  may  be  formulated  as 
an  optimization  problem  with  the  objective  function 


Z 


(1-5) 


J  o 

if  maY  be  measured  as  a  continuous  function  of  time.  If  the 

A 

measurements  y_  are  made  at  discrete  values  of  the  independent 
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variable,  then  the  objective  function  is 

2 

Z  I  | | V i  -  V ±  ||  (1-6) 

i 

With  the  first  objective  function  (1-5)  the  identification 
problem  has  been  formulated  as  an  optimal  control  problem. 
When  the  second  objective  function  must  be  used,  the 
identification  problem  has  been  formulated  as  a  multi¬ 
point  boundary  value  problem.  It  is  this  formulation  and 
the  solution  of  the  resultant  boundary  value  problem  that 
is  the  subject  of  the  remainder  of  this  work. 


QO 
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II .  THEORY 

A .  Formulation  of  the  Problem 

Suppose  that  the  behavior  of  a  physical  system  can 
be  adequately  approximated  by  the  following  vector  differen¬ 
tial  eauation: 


dx 

—  =  f  (x,  £,  t)  (II-l) 

dt 

where  x  =  vector  of  state  variables  (order  p) 

a  =  vector  of  constant  system  parameters  (order  q) 
t  =  independent  variable 


Suppose  further  that  n  different  experiments  have 
been  conducted  on  the  actual  physical  system  with  the  initial 
conditions  x_.  (0)  ,  j  =  1,  2,  ...  n.  At  certain  discrete 

values  of  t,  some  or  all  of  the  elements  of  x  have  been 
observed  and  recorded  as  data. 

The  set  of  vector  differential  equations  describing 
all  these  experiments  is: 

dx  . 

— 1  =  f  •  (x.,  a ,  t) ,  j  =  1,  2,  ...  n  (II-2) 

dt  D  D 


with  the  initial  conditions  x_.  (0)  and  the  functions 

fj  (x^,a,  t)  being  specified  for  each  experiment.  The  data 

pertinent  to  the  j  experiment  will  be  referred  to  as  the 
+• 

j  data  set  throughout  the  remainder  of  this  work. 
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th 

If  is  taken  to  be  the  k  discrete  value  of 

t  for  the  j  experiment,  then  x^ .  ^  may  be  defined  as  the 


.  th 


observed  value  of  the  i  state  variable  at  t.,  . 

Dk 

Since  x^  represents  the  value  predicted  by  equation 

(II-2)  for  the  value  t.,  of  the  independent  variable,  an 

error  e.  may  be  defined  as: 

13K 


e  .  . ,  =  x..,  -  x... 

ljk  ijk  ijk 


(II-3) 


Let  £  be  defined  as  the  error  vector  containing  all  the  ele¬ 
ments  ordered  in  some  arbitrary  way. 

The  problem  may  now  be  posed  mathematically  as  one 
of  finding  a  such  that  some  objective  function  Z  of  £  is 
minimized.  The  two  most  popular  objective  functions  are 
1)  the  "least  squares"  criterion: 


(II-4) 


where  W  is  a  positive  definite  diagonal  matrix  with  the 
diagonal  elements  representing  weighting  factors,  and  2) 
the  "Chebyshev"  criterion: 


Z 


max 


ijk  ,Wijk  °ijk 


(II-5) 


To  solve  this  problem  directly  as  it  stands,  an 
optimum  seeking  method  would  be  required  such  as  Rosenbrock's 
hill  climbing  procedure  (25)  or  the  Marquardt  (22)  procedure 
as  proposed  by  Ball  and  Groenweghe  (1) .  Another  way  to  solve 
this  problem,  as  shown  by  Bellman  and  Kalaba  (2),  is  to 


i  +<§ 
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formulate  it  as  a  boundary  value  problem  by  taking  the  para¬ 
meter  vector  a  to  be  the  solution  of  the  vector  differential 
equation : 

da 

—  =0  ( I I- 6 ) 

dt 


Now  the  complete  system  of  differential  ecmations  may  be 
written  as: 


dx  . 
-3 

dt 


da 

dt 


f j  (x.,  a,  t)  j  =  1,  2  ...  n 


(II-7) 


0 


The  initial  conditions  x_.  (0)  are  specified  for  j  =  1, 

2  ...  n  but  a  (0)  is  unknown. 

This  is  inherently  a  boundary  value  problem,  as  the 
initial  condition  for  a  (0)  is  not  known.  Equation  (II-7) 
may  be  rewritten  in  a  more  condensed  form  by  letting: 


Z  = 


and 


g  (y /  t) 


fq  (xq/  a,  t) 

f  .  (x  ,  a,  t) 
— n  — n  — 

0 


(II-8) 
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Then,  in  the  new  notation,  the  problem  becomes: 

—  =  2.  (£,  t)  ( XI-9 ) 

dt 

Y_  and  2.  are  vectors  of  order  r  x  1  where  r  =  np  +  a. 

The  problem  of  determining  the  optimum  parameters 
which  characterize  the  given  physical  system  has  now  been 
formulated  as  a  boundary  value  problem.  This  boundary 
value  problem  is  almost  certainly  nonlinear. 

B .  Quasilinearization  and  the  Solution  of  Nonlinear  Boundary 
Value  Problems 

There  are  no  good  general  techniques  for  handling 
nonlinear  boundary  value  problems.  However,  there  exists 
powerful  techniques  for  linear  boundary  value  problems  which 
make  their  solution  almost  routine.  Bellman  and  Kalaba  (2) 
have  taken  advantage  of  this,  by  using  Kantorovich's  (13) 
extension  of  the  Newton-Raphson  (6)  method  to  function  space, 
which  results  in  a  succession  of  linear  boundary  value  prob¬ 
lems  which  hopefully  converge  to  the  original  nonlinear 
boundary  value  problem. 

The  seauence  of  linear  problems  is  defined  by  the 
following  recurrence  relationship: 
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dz 


(m+1) 


dt 


=  2(m)  +  jm 


(m+l)  (m) 

£  -  1 


=  j<m>  y<m+1>  + 


_a(m)  -  J(m>  y<m>  ] 


m  =  0 ,  1 ,  2 ,  ... 


(11-10) 


where  J  ^  (t) 


is  the  Jacobian  matrix  of  partial  derivatives. 


(m) 


3gi  _ 

9  g 

•  •  • 

3  Yl 

3  y 

9  g 
^r 

3g 

•  •  • 

3  Y1 

3  y 

(11-11) 


Given  an  initial  guess  solution  and  solving  the  succes¬ 

sion  of  linear  problems  represented  by  equation  (11-10)  with 
the  same  boundary  conditions  as  given  for  the  original  problem, 
the  recurrence  relations  will  converge  to  the  solution  of  the 
nonlinear  problem  such  that: 


limit 

m  -*  «> 


(H-12) 


To  assure  convergence  of  this  scheme  for  any  arbitrary 
initial  guess  ,  certain  limitations  on  equation  (II-9) 

must  be  imposed.  First,  each  of  the  functions  which  make  up 
the  function  vector  g(y_,  t)  must  be  strictly  convex.  By 
definition,  this  means  that  for  each  function,  the  Hessian 
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matrix  is  positive  definite.  The  Hessian  matrix  is  defined 
as : 


.  2 
3  g 


3  x  .  3  x  . 
_  i  D 


3  2g} 
3  2x. 


3  2g, 


3  X,  3  X 
1  r 


k  =  1 ,  2  .  .  .  r 


3\ 

3  x,3  x 
1  r 


3  2g, 


3  2x 


r  -1 


(H-13) 


This  property  makes  it  possible  to  write: 

dZ 

—  =  £  (v,  t)  =  Max  [_g(z,  t)  +  J  (z ,  t)  (y-^Q  (H-14) 

dt  z 


Now  if  the  function  w  is  introduced  as  the  solution  to  the 
associated  linear  equation: 


dw 

—  =  £  (z,  t)  +  J  (z^,  t)  (w-z)  (II-]4) 

dt 


then,  provided  that  a  certain  positivity  property  holds,  it 
has  been  shown  by  Kalaba  (12)  that  the  solution  of  (11-14) 
is  given  by: 
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Y_  =  Max  w  (t,  z)*  (11-15) 

z 

Once  this  has  been  established,  Kalaba  (12)  has  further 
proved  that  the  sequence  of  linear  problems  represented  by 
equation  (11-10)  will  converge  monotonically  and  quadra- 
tically  to  the  solution  of  the  original  nonlinear  problem. 
The  required  positivity  property  is  that  if 


du 

—  -  J  (z,  t)  u  =  6  >  0 
dt  _  _  _  _ 


with  u  (0)  =  0  (11-16) 

then  u  >_  0,  t  >_  0 

A  sufficient  condition  for  this  property  to  hold  is  that  the 
off-diagonal  elements  of  J  should  be  non-negative  as  shown 
by  Kalaba  (12) . 

The  convexity  and  positivity  properties  discussed 
above  are  not  present  in  a  large  number  of  practical  problems . 
However,  this  algorithm  will  often  converge  in  spite  of  this, 
provided  that  the  initial  guess  is  sufficiently  close  to 

the  answer 


* 


A  similar  minimization  operation  can  be  defined  for 
concave  problems 


' 
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To  obtain  an  initial  guess  within  the  domain  of 
convergence  of  a  particular  problem  is  often  no  trivial 
task.  An  obvious  approach  would  be  to  use  some  type  of 
multi-variable  search  or  hill-climbing  method  to  obtain  an 
approximate  result. 

Another  less  obvious  approach  would  be  to  set  up  an 
artificial  problem  whose  solution  is  closer  to  the  arbitrary 
initial  guess.  This  is  the  approach  followed  in  the  remainder 
of  this  work. 


C.  Boundary  Value  or  Data  Perturbation 

Let  the  initial  guess  be  the  solution  of  the 

following  vector  differential  equation: 


(0) 


=  2  (y 


(O) 


t) 


(II-17) 


dt 

with  the  initial  condition  y^  (0)  =  b 

where  the  constant  vector  b  is  an  intelligent  guess  of  the 

actual  initial  condition  y  (0)  which  is  not  entirely  known. 

„  th 

Now  let  y^  be  the  given  boundary  value  for  the  i 

element  of  y  at  a  discrete  value  of  the  independent  variable 

t^  and  let  y!^  be  the  value  of  y^k  predicted  by  equation 

(16)  . 


If  is  within  the  domain  of  convergence  for  the 

given  problem,  the  quasilinearization  procedure  may  be  applied 


•  ' 


- 


.. 

..  • 

:  :  '  ,  /  n  1 
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directly.  However,  if  the  initial  guess  y^^  is  outside  the 
domain  of  convergence  for  the  actual  problem,  it  should  be 
in  the  domain  of  convergence  for  another  problem. 

If  an  artificial  set  of  boundary  conditions  y£^  are 
defined  by: 


ik 


y  (0) 
yik 


R(yik  -  vi?*) 


ik 


0  <  R  <  1 


(11-18) 


then  perhaps  the  quasilinearization  procedure  may  be  applied 
to  solve  the  derived  problem.  Once  the  derived  problem  is 


solved,  the  solution  y*  may  be  used  as  the  initial  guess  £ 


(0) 


and  the  procedure  may  be  repeated  until  y^0^  is  sufficiently 
close  to  the  actual  solution  so  that  the  actual  boundary 
values  y,^  may  be  used  in  place  of  the  pseudo-boundary  condi¬ 
tions  y+^.  This  stepping,  as  illustrated  in  Figure  1  should 
always  succeed,  provided  there  is  a  domain  of  convergence 
around  y/°^  which  is  sufficiently  large  so  that  the  solution 
may  be  reached  in  a  finite  number  of  steps. 

Applying  this  procedure  to  the  parameter  identification 
problem  is  straight  forward.  The  only  initial  conditions  that 
are  not  known  are  those  for  the  parameter  vector.  Therefore, 
to  define  y/°^ ,  it  is  necessary  only  to  supply  a  guess  for 
each  of  the  unknown  parameters.  The  data  is  treated  just  as 
boundary  values  in  the  perturbation  stepping  procedure. 


' 
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D .  The  Linear  Boundary  Value  Problem 

A  systematic  method  for  solving  the  successive 
linear  boundary  value  problems  will  now  be  discussed. 

Each  of  the  linear  boundary  value  problems  is  of 

the  form: 


dy 

dt 


(m+ 1 ) 


,  (m)  (m+1)  ,  (m) 

=  J  y  +  g 


y(m)  (11-18) 


with  the  appropriate  boundary  conditions.  The  vectors  g 


(m) 


/  \  / wn  \ 

and  vl  ,  and  the  matrix  J  '  are  all  known  functions  of  the 
independent  variable  t.  The  general  solution  of  (11-18)  is: 


^(m+1)  =  y (m+1)  c (m+1 )  + 


u 


(m+1 ) 


(11-19) 


( m+ 1 ) 

as  given  in  Coddington  and  Levinson  (4) .  X(t)  is  the 
fundamental  matrix  which  is  the  solution  to  the  matrix  differ¬ 


ential  equation: 


dY 


(m+1) 


_  j (m)  y (m+1) 


(11-20) 


dt 


/  ry.  -I  \ 

with  the  initial  conditions  =  l^rxrj 

The  vector  uj™+^  is  a  particular  solution  to  the  vector 
differential  equation: 


du 


(m+1) 


dt 


=  j(m)  H<m+1)  +  a(m)  _  j(m)  ^(ra) 


(11-21) 


. 
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with  the  initial  condition 


The  constant  vector  c 

(m+1 ) 


(m+1) 


is  the  vector  of  the  initial 


condition  of  ^(0)  which  may  be  determined  from  the  boundary 
conditions . 

In  principle,  the  above  relations  may  be  applied 
directly  to  any  boundary  value  problem  using  a  numerical 


v (m+1) 
-(t) 


and  u  •  However, 


integration  procedure  to  generate  Y 


if  a  large  number  of  experiments  are  to  be  considered,  each 

with  a  different  initial  condition,  and  if  the  number  of 

state  variables  is  large,  then  the  order  r  of  the  vector  £, 

as  can  be  seen  from  equation  ( 1 1— 8 ) ,  is  quite  high.  The  total 

number  of  differential  equations  that  would  have  to  be  inte- 

2 

grated  numerically  is  r  +  r.  The  value  of  r  is  np  +  q  where 
n  is  the  number  of  experiments,  p  is  the  number  of  state 
variables,  and  q  is  the  number  of  parameters. 


However,  J ^  has  a  special  structure,  and  because 


of  this,  the  magnitude  of  the  problem  may  be  reduced  consider- 


/  m  \ 

ably.  The  structure  of  J  ,  as  can  be  seen  by  again  referring 
to  equation  (II-8),  is  as  follows: 


ViBbnuod  erid  mD'it  beciltmo^&b  sd  vsm  do idw  v  i 
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0 


(m) 


with 


j(m) 

~XD 


and 


0 


0 


J  . 
-a: 


O 

0 


J  J 
— xn  —an 

O  0 


3 


3  x  . 
P3 


3  f 


!£2 


*lj 


3  f  . 

R1 


3  x 


p: 


j  =  1 ,  2 ,  ...  n 


3  f 


li 


3  a. 


3  f  . 

,  .  PI 

3  a. 


qxi 

3  a 


3  f  . 

El 

3  cl 


(m) 


(11-22) 


(11-23) 


(11-24) 


1,  2, 


•  •  • 


n 
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The  equation  (11-18)  may  be  rewritten  as: 


d 


dt 


1 

1 — 1 

xl  • 

! 

(m+1) 

1 

•  I* 

i— 1 

_ 1 

(m+1) 

1 

•  |  i-h 

1— 1 

_ 1 

(m) 

r 

• — 1 

xl  • 

1 _ 

• 

• 

X  . 

=  j(m) 

• 

• 

x  . 

~3 

• 

• 

+ 

• 

• 

f  . 

~3 

• 

• 

-  jW 

• 

X  . 

~3 

• 

• 

~3 

• 

• 

• 

X 

• 

X 

• 

f 

• 

X 

— n 

— n 

—n 

— n 

a 

a 

0 

a 

(m) 


(11-25) 


/  m4_  "I  \ 

Since  each  vector  x.  'is  coupled  only  with  the  parameter 
vector  a  t  and  not  with  any  other  state  variable  vector 

from  a  different  experiment,  the  general  solution  can  be 
written  as: 

(m+1) 


si 

(m+1) 

F  ,  . 

— xl 

•  o 

Sal 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

•  • 

• 

• 

x  . 

— 

0  F  . 

0 

G  . 

~3 

~x3 

— 

-a: 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

X 

0 

F 

G 

— n 

a 

— xn 

0 

—an 

F 

—a  _ 

(m+1) 

r-  — 

(m+1 ; 

-1 

• 

“-1  ~ 

• 

• 

• 

c  . 

+ 

• 

u  . 

-3 

• 

• 

-3 

• 

• 

• 

c 

• 

u 

— n 

— n 

c 

u 

a_ 

_  — a  _ 

(11-26) 


where  the  first  term  on  the  right  hand  side  represents  the 
fundamental  matrix  and  the  last  term  represents  a  particular 


1 


. 

■ 
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solution . 


The  block  diagonal  entries  in  the  fundamental  matrix 


F  (m+1) 

-xj 

(3  1/  2 ,  ...  n) 

and  F(m+1) 
—a 

are  the  solutions 

of ; 

dF(m+l)  _  (m) 

-xj  -xj 

F<m+1)  (j 

— X  j  J 

=  1 ,  2 ,  ...  n) 

dt 

p (m+i )  (o)  =  j 

-xj  -(pxp) 

(11-27) 

V 

and 

dF(m+i) 

^  -0' 

P (m+1) 
—a 

(0)  =  I  ,  . 

— ( qxq ) 

(11-28) 

dt 


From  (11-28)  it  follows  that 


p  (m+i)  (  j  _ 

-a  -(qxq) 


The  other  block  entries  in  the  fundamental  matrix  G 
(j  =  1,  2  ...  n)  are  the  solutions  of: 


(11-29) 

(m+1) 

-aj 


dG(m+1) 

_ 

dt 


j  (m)  G  (m+1 )  +  j  (m)  F(m+1) 

— xj  — aj  — aj  —a 


j (m)  (m+1)  +  j (m) 

-xj  -aj  -aj 


(11-30) 


with  the  initial  condition 


(m+i)  (0)  =  0 

-aj  -(pxq) 


The  entries  in  the  particular  solution  u^  (j  =  1,  2.  ...  n) 

and  u  are  the  solutions  of : 

—a 
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du!m+1) 

~3 

dt 


J(m) 

-X} 


u(m+l)  +  (m) 

~3  -a  3 


u(ra+1)  +  f. 
-a  -j 


J(m) 

~XD 


(11-31) 


with  the  initial  condition  ufm+^^  (0) 


~3 


-(pxl) 


and 


du  <m+1) 
—a 

dt 


=  0 


(11-32) 


with  the  initial  condition  (0)  =0, 

-(qxl) 


—a 


(11-33) 


and,  therefore ,  equation  (11-31)  may  be  rewritten  as: 


du 


dt 


(m+1) 

2  =  j(m)  u(m+l)  +  f  _  x(m)  _  j  (™)  a(m) 


~XD  ~3 


— x 


with  the  initial  condition  u^m+(0^)  =0,  ,  x. 

~3  -(pxl) 


-a:  - 

(11-34) 


Thus,  the  number  of  differential  equations  to  be 

2 

integrated  has  been  reduced  to  (np  +  npq  +  np)  from 
2 

(np  +  q)  +  (np  +  q) .  This  is  a  considerable  reduction  in 
the  memory  and  computation  requirements  for  this  method. 


Equation  (11-26)  may  be  expanded  in  terms  of  each 

(] 
j 


state  variable  vectors  xfm+^  as  follows: 


x  (m+1 )  _  F  (m+1)  c  (m+1 )  +  G(m+1)  c(m+l)  +  U(™+1) 

-j  -xj  -j  -aj  -a  -j 

(11-35) 


It  is  known  from  the  formulation  of  the  problem  that: 


. 


. 
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c!m+1> 

=  x!m+1) 

(0) 

~3 

-3 

c(m+l) 

—a 

=  a (m+1) 

(0) 

(11-36) 

=  a(m+1) 

In  a  large  number  of  problems,  Xj  (0)  is  known,  and  there¬ 
fore  can  be  assumed  to  be  exact.  This  assumption  further 
simplifies  the  problem. 


Let 


(m+1 ) 
v . 

“3 


be  defined  by: 


v  (m+1 )  _  p  (m+1)  c(m+l)  +  u  (m+1) 

~j  -xj  -j  -j 


Using  this  definition,  equation  (11-35)  becomes: 
x  (m+1)  _  v  (m+1 )  +  Q  (m+1 )  a(m+l) 

-j  ~j 


(H-37) 


(11-38) 


With  the  assumption  that  x_.  (0)  is  exact,  the  integration 

of  equation  (11-27)  and  (11-31)  may  be  replaced  with 

v !m+1^  which  is  the  solution  to: 

-3 


dv) 

=1 


(m+1) 


dt 


j(m) 

-xj 


(m+1)  (m) 

v  '  -  x  ; 

~3  ~3 


(m) 


J(m)  a <m)  +  f  ; 
~a3  -  -3 


(11-39) 


with  the  initial  condition  (0)  =  x_.  (0)  (j  =  1,  2 


n) 


The  matrix  must  still  be  evaluated  as  before, 

—a  3 

with  equation  (11-30).  Thus,  all  the  vectors  (j  =  1, 

2  ...  n)  can  be  expressed  linearly  in  terms  of  a  (m+-^  alone 


- 
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and  the  least  squares  analysis  may  be  carried  out. 

This  simplification  means  that  the  number  of  differ¬ 
ential  equations  per  iteration  is  only  (np  +  npq)  as  com- 
2 

pared  to  (np  +  npq  +  np)  for  the  case  where  the  initial 
conditions  are  treated  as  unknown.  The  number  of  variables 
to  be  considered  in  the  least  squares  analysis  has  also  been 
decreased  from  (np  +  q)  to  (q) .  This  is  impofctant  because 
of  the  inherent  ill-conditioning  problems  associated  with  the 
least  squares  analysis.  This  problem  can  also  be  avoided  by 
using  the  Chebyshev  criteria. 
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III.  PRACTICAL  APPLICATION 


The  theoretical  concepts  presented  in  the  previous 
section  can  be  incorporated  into  a  practical  algorithm. 

To  do  this,  a  step  by  step  procedure  must  be  outlined  in 
an  organized  way,  which  may  be  programmed  for  the  digital 
computer . 


A.  Data  Perturbation  and  a  Stepping  Procedure 


An  empirical  procedure  was  devised  for  the  purpose 
of  controlling  the  step  size.  Perhaps  a  more  efficient 
method  could  be  devised.  However,  this  one  has  been  success¬ 
fully  implemented.  The  procedure  is  outlined  below. 

Step  1:  The  maximum  step  size  S  ,  the  maximum  allowable 

c  max 

relative  change  in  the  parameters  for  one  iteration 

C  ,  and  the  minimum  allowable  change  in  the  para- 
max  '  ^ 

meters  for  one  step  C  .  are  all  specified.  The 

min 


Step  2: 


Step  3: 


initial  guess  for  the  parameter  vector  a 
supplied . 


(0) 


is 


,max 


The  maximum  element  in  the  error  vector  ( i i &i , j ,k 
is  computed  using  the  current  value  of  a ^ . 


If  max 


i,j,k  ijk 


>  S 


max 


R  is  defined  by  the  equation 

S 


max 


R  = 


max 

i, j ,k  i, j  ,k 


1 
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Step  4: 


Step  5: 


If  max 
i,  j  fk 


e  .  .  . 

if  D  /k 


<  S 

-  max 


then  R  =  0 

Define  the  pseudo-data  set  as 


X .  .  ,=  X .  .  .  +  R  (X.  .  , 

if Dfk  if ] ,k  if]fk 


x 


if  j  fk 


) 


for  all  values  of  i,  j,  k 

Note  that  if  R  =  0  then  the  pseudo-data  set  is 
equal  to  the  real  data  set. 

The  quasilinearization  procedure  is  initiated  with 

a^°)  being  the  initial  guess  for  the  parameter 

* 

vector  and  the  data  set  x.  .  ,  (for  all  i,  j,  and 

1  t  3  r  ■K 

k)  being  the  required  boundary  conditions.  Norm¬ 
ally,  this  STEP  is  terminated  if  some  convergence 
criteria  is  met  or  if  a  specified  number  of  itera¬ 
tions  has  been  completed.  If  during  the  procedure 


q 

l 

i=l 


m+1 
a . 
i 


m+1 
a . 

l 


>  C 

max 


for  any  m  (m=l  /  2,  ...)  then  the  procedure  is 

terminated,  S  is  halved  and  STEP  2  is  re-initiated. 

max 

If  R  =  0  and  STEP  5  is  successfully  completed,  the 
problem  is  solved  and  the  perturbation  procedure 
may  be  terminated. 


Step  6 : 


■ 
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Step  7: 


Step  8: 


If 


q 

l 

i=l 


.  (M) 

d  . 

1 


(0) 


/ 


a.(M) 

l 


<  C 


mift 


then  S  is  tripled, 
max  ^ 

a/0^  is  defined  as  being  equal  to  a^  and 
Step  2  is  re-initiated. 


The  above  outline  covers  the  basic  steps  used  in 
attempts  to  solve  the  indentif ications  problem.  Many  of 
the  finer  details  are  not  included  here  for  the  sake  of 
simplicity.  The  complete  algorithm  is  listed  in  the  form 
of  a  Fortram  IV  program  in  Appendix  A. 


B .  Quasi linearization 

In  principle,  the  data  perturbation  procedure  does 
not  require  the  use  of  quasilinearization.  It  does,  how¬ 
ever,  require  a  rapidly  converging  method  for  the  solution 
of  the  nonlinear  boundary  value  problems  encountered  with 
each  step.  As  the  perturbation  insures  that  the  initial 
guess  is  quite  close  to  the  solution  for  each  step,  quasi¬ 
linearization  is  ideally  suited. 

Assuming  the  initial  conditions  to  be  exact  is  not 
essential  to  data  perturbation.  This  assumption  is  made 
solely  to  reduce  the  computational  load  in  the  evaluation 
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of  the  fundamental  matrix  and  in  the  subsequent  least  squares 
analysis . 

To  solve  each  of  the  successive  nonlinear  boundary 
value  problems,  the  following  equations  must  be  integrated 
simultaneously  for  each  iteration. 


dx 

-1 


(0) 


dt 


3  1,  2,  ...  n 


(III-l) 


dt 


(k-1 )  ~1 
x  . 

-1 


+  f  .  (x : 
~D 


(k) 


t) 


j  =  1,  2,  . . .  n  (III-2) 

k  ~~  1,  2,  ...  m 


dv!ra+1) 

_ 

dt 


(m)  (m+1) 

U  V  . 

-XI  -1 


J(m) 

-a: 


a 


(m) 


+  f(m+l) 

-j 

(x(ra),  a(m,)  , 

j 

dG(m+1)  ,  . 

=  1 ,  2 ,  ...  n 

c  (m+1)  +  j  (m) 

dt  -X^> 

-aj  -aj 

j 

=  1 ,  2 ,  ...  n 

(HI-3) 


( 1 1 1  —  4  ) 


initial  conditions 


(0) 
x  . 

-1 

(0)  =  x. 

(0) 

j  =  If 

2, 

...  n 

(k) 
x  . 

—J 

(0)  =  x_. 

(0) 

j  = 

k  =  1, 

2, 

2, 

...  n 

.  .  .  m 

with  the 


* 
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v]m+1)  (0)  =  Xj  (0)  j  =  1,  2,  ...  n 

G(m+1)  (0)  =  o  j  =  1,  2,  ...  n 

J 

Equation  (III-l)  is  used  to  generate  the  initial 
guess;  equation  (III-2)  is  used  to  generate  the  previous 
m  solutions;  and  equations  (III-3)  and  (III-4)  are  required 
to  determine  . 

All  the  parameter  vectors  a^0/ ,  a  ^  ,  ...  am  must 

be  stored  to  be  used  in  the  evaluation  of  the  above  system 
of  differential  equations. 

Once  the  above  system  is  solved,  the  following 
relationship 

G(r!l+1)  a  +  vfm+1)  -  xfm+1)  =  0  (III-5) 

— a  j  -  -3  -3 

may  be  used  to  form  the  linear  constraints  required  for 
the  solution  of  the  least  squares  problem.  In  general 
these  constraints  may  not  always  be  linear. 

To  extend  the  algorithm  to  include  nonlinear 
boundary  conditions  would  not  be  difficult  in  principle. 
However,  it  was  decided  to  leave  the  consideration  of  more 
complex  boundary  conditions  to  be  studied  in  future  work. 


■ 


' 


; 
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IV.  APPLICATION  IN  CHEMICAL  REACTION  KINETICS 

At  the  heart  of  almost  any  major  chemical  process 
is  the  reactor.  Reactor  engineering  and  design  is  thus  of  un¬ 
questionable  importance  to  the  chemical  industry. 

To  design  a  reactor  effectively  requires  an  intimate 
knowledge  of  the  chemical  kinetics,  and  it  is  highly  desirable 
to  have  some  type  of  reliable  mechanistic  rate  equation.  The 
determination  of  kinetic  rate  parameters  has,  justifiably,  been 
of  prime  concern  in  the  analysis  of  kinetic  data. 

E.S.  Lee  (18,  19)  has  suggested  the  application  of 
quasilinearization  to  the  identification  of  kinetic  rate  para¬ 
meters.  He  has  illustrated  this  approach  by  considering  the 
difficult  problem  of  analyzing  the  kinetic  data  obtained  from 
a  tubular  reactor  in  which  axial  diffusion  was  important.  He 
considered  both  nonlinear  boundary  conditions  and  the  analysis 
of  non-isothermal  data.  Unfortunately,  however,  quasiiineari-- 
zation  is  unstable  for  many  problems. 

If  this  is  the  case  for  a  particular  problem,  either 
the  problem  must  be  modified,  or  some  way  must  be  found  to 
extend  the  domain  of  convergence.  Boundary  value  or  data 
perturbation  is  one  hope  fordoing  the  latter. 

A .  A  Nonlinear  Example 

A  hypothetical  reaction  described  by  the  following 
differential  equations  was  studied. 


- 


30 


dx 


1 


dt 


-  a1  (x1  -  kelx2)  -  a2 (x1  -  kg2  x3) 


dx. 


dt 


=  al (xi  -  keix2>  "  a3(x2  -  ke3x3> 


(IV-1) 


Following  the  conventions  set  forth  earlier,  the 
state  vector  is 


( IV- 2 ) 


with  x^  and  x2  being  the  composition  of  components  one  and 
two.  The  composition  of  component  three  is  given  by  the 
following  linear  relationship: 

x^  =  1  -  x1  -  x2  ( IV- 3 ) 


The  unknown  parameter  vector  is 


( IV- 4 ) 


The  other  parameters  k  ^ ,  ke2,  and  ke2  are  equilibrium  con¬ 
stants  which  are  assumed  to  be  known. 

The  purpose  for  considering  this  example  was  to  test 
out  the  algorithm  and  to  examine  the  convergence  rates.  To 
do  this,  values  were  chosen  for  the  parameter  vector  and  the 
equilibrium  constants  and  exact  data  was  generated  for  three 
different  initial  conditions  as  given  in  Table  1.  The  parameter 
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and  equilibrium  constants  used  to  generate  this  data  were: 


a^  —  2.0 

k  .  =  1.8 
el 

a2  =  3,5 

k  0  =  3.0 
e2 

a3  =  5.0 

o 

• 

i — 1 

ii 

co 

<U 

Using  the  algorithm  developed  in  the  previous  section,  the 
parameter  vector  was  then  determined  from  the  data  in  Table 
1. 

In  order  to  evaluate  the  equations  (III-2) ,  (III-3) , 

and  (III-4),  the  following  relations  are  used  together  with 
an  initial  guess. 


J  . 
-*D 


J  . 
-a  3 


2a1x1  +  a2  (1  +  ke2) 


2alxl 


a-jk  _ 

3  e3 


a^k 


a^k 


el 


el 


j  =  1,  2,  3 


2 

Ui 

-  k  ,  x~ 
el  2 

- 

X1  +  ke2  (1 

_  2 
X1 

-  k  ,x_ 
el  2 

L°I 

-  x0  +  k  - 
2  ei 

+  a3  (1  +  ke3) 

( IV- 6 ) 


(1  -  xx  -  x2) 


j  =  1,  2,  3 


(IV-7) 


. 

Data  for  the  Nonlinear  Example 
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f  . 
-3 


kelx2> 


k  -.x-, 
e3  3 


kelx2> 


ke3x3> 


( IV- 8 ) 


j  =  1,  2,  3 


xx  (0) 


x2(°) 


x3(0)  = 


(IV- 9 ) 


The  first  guess  vector  considered  was 


0.0000001 

0.0000001 

0.0000001 


(IV-10) 


From  this  initial  guess,  the  quasilinearization  procedure 
could  be  used  directly  without  the  aid  of  data  perturba¬ 
tion.  The  successive  parameter  vectors  are  given  in  Table 
2.  The  convergence  of  the  corresponding  concentration  pro¬ 
files  is  shown  in  Figures  2,  3,  and  4.  The  rapid  conver¬ 
gence  of  the  quasilinearization  procedure  is  verified. 
However,  as  mentioned  before,  the  domain  of  convergence  is 
often  quite  small. 
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TABLE  2 


Convergence  of  the  Nonlinear  Example 
Without  Data  Perturbation 


Iteration 

ai 

a2 

a3 

0 

0.00 

0.00 

0.00 

1 

0.48 

0.36 

0.51 

2 

1.16 

1.01 

1.68 

3 

2.01 

1.92 

3.08 

4 

2.26 

2.75 

4.13 

5 

2.06 

3.33 

4.83 

6 

2.00 

3.49 

5.00 

CONCENTRATION 
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TIME 

Figure  2  Convergence  of  nonlinear  example  using 
Quasilinearization  (data  set  one) 


CONCENTRATION 
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TIME 

Figure  3  Convergence  of  nonlinear  example  using 
Quasilinearization  (data  set  two) 


CONCENTRATION 
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TIME 

Figure  4  Convergence  of  nonlinear  example  using 
Quasilinearization  (data  set  three) 
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With  another  initial  guess. 


“  10“ 
10 
10 


(IV-11) 


attempts  to  find  the  solution  using  the  quasilinearization 
directly  were  frustrated.  In  one  iteration,  the  parameter 
vector  changed  from  the  values  given  in  (IV-11)  to 


a 


(1) 


“-2.74“ 

-73.88 

-17.12 


(IV-12) 


With  the  values  in  (IV-12)  the  integration  of  the  funda¬ 
mental  matrix  becomes  unstable  and  further  progress  was 
not  possible. 

However,  by  perturbating  the  data  to  within  ten 
percent  of  the  concentration  profiled  generated  from  the 
initial  guess,  a  solution  can  be  obtained  through  the 
stepping  procedure  combined  with  quasilinearization.  After 
five  steps,  the  derived  initial  guess  is  within  the  domain 
of  convergence  for  the  direct  application  of  the  quasilineari¬ 
zation  procedure.  The  stepping  procedure  produced  the  sequence 
of  derived  initial  guess  given  in  Table  3.  For  each  step, 
a  maximum  of  three  quasilinearization  iterations  were  applied 
to  determine  the  next  initial  guess.  If  convergence  to 
three  significant  figures  was  obtained,  the  step  was  consider¬ 
ed  completed. 
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The  successive  profiles  yielded  by  the  stepping 
procedure  are  shown  in  Figures  5,  6,  and  7. 


TABLE  3 

Convergence  of  the  Nonlinear  Example  with 
Data  Perturbation 


Step  No. 

ai 

a2 

a3 

0 

10 

10 

10 

1 

9.61 

7.96 

9.32 

2 

8.36 

5.64 

8.27 

3 

7.28 

4.76 

7.65 

4 

4.83 

3.83 

6.40 

5 

3.12 

3.58 

5.54 

' 
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1.05 


0.00 


0.00 


0.02 


0.04  0.06 

TIME 


0.08 


0.10 


Figure  5  Convergence  of  nonlinear  example  using 
Data  Perturbation  (data  set  one) 


CONCENTRATION 
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Figure  6  Convergence  of  nonlinear  example  using 
Data  Perturbation  (data  set  two) 


CONCENTRATION 
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TIME 

Figure  7  Convergence  of  nonlinear  example  using 
Data  Perturbation  (data  set  three) 
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With  the  completion  of  the  stepping  procedure,  the 
initial  guess  is 


a 


(0) 


3.12 

3.58 

5.54 


( IV- 13 ) 


From  this  guess,  the  direct  application  of  the  quasilineari¬ 
zation  procedure  yielded  the  answer  auite  rapidly  as  shown 
in  Table  4 . 


TABLE  4 

Convergence  of  the  Last  Step 
for  the  Nonlinear  Example 


Iteration  No. 

ai 

a 

2 

a3 

0 

3.12 

3.58 

5.54 

1 

1.81 

3.45 

4.89 

2 

2.002 

3.498 

5.000 

3 

2.000 

3.500 

5.000 
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From  this,  it  is  evident  that  data  perturbation, 
together  with  quasilinearization,  can  be  a  powerful  tool 
in  the  solution  of  this  type  of  problem.  The  main  dis¬ 
advantage,  at  present,  is  that  there  is  no  way  that  the 
optimum  step  size  can  be  determined  in  advance.  Choosing 
it  too  small  will  waste  a  large  amount  of  computational  time, 
but  it  must  be  chosen  small  enough  to  ensure  convergence. 

At  present,  the  empirical  procedure  for  changing  the  step 
size  mentioned  in  the  previous  chapter,  is  used.  This  pro¬ 
cedure  has  proven  satisfactory  in  this  problem  as  well  as  the 
others  in  this  work.  It  is  recommended,  but  there  is  room 
for  a  considerable  amount  of  improvement. 

B .  Sensitivity  to  the  Choice  of  Experiments  and  to  Experimental 
Error 


Wei  and  Prater  (27)  presented  a  method  for  identify¬ 
ing  the  kinetic  rate  constants  in  a  linear  kinetic  model. 

This  method  involved  choosing  experiments  in  an  organized 
manner  so  that  the  eigenvector  directions  characteristic  of 
the  system  could  be  determined.  Once  all  the  eigenvector 
directions  are  determined,  the  rate  constants  may  be  calcula¬ 
ted  directly.  If  some  of  the  eigenvector  directions  are  not 
known,  the  rate  constants  can  not  be  uniquely  determined. 
Logically,  it  follows  that  the  experiments,  from  which  the 
data  sets  are  taken,  must  have  strong  components  in  each  of 


. 

. 
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the  eigenvector  directions  if  the  parameter  vector  is  to  be 
uniquely  determined.  It  would  also  seem  probable  that  an 
analogous  situation  would  exist  with  nonlinear  systems. 

To  test  these  principles,  an  example,  which  was  used 
be  Wei  and  Prater  (27) ,  was  chosen. 


1  -  butene 


The  above  mechanistic  diagram  may  be  represented 
by  the  following  system  of  differential  equations. 


dt 

dx2 

dt 


_al  <X1  '  k31x2>  '  a2  (X1  '  ke2x3) 


( IV- 14 ) 


kelx2) 


a3  (x2  “  ke3X3) 


with 

cis 

2  - 


x^  and  x^  being  the  concentration  of  1-butene  and 
-  2  -  butene  respectively.  The  concentration  of  trans  - 
butene  which  is  x^  is  given  by  the  relationship 


X3 


x. 


1 


xl“ 


( IV- 15 ) 


■ 
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The  equilibrium  constants  are: 


k 


el 


k 


e3 


0.4469 


0.2685 


0.6002 


(IV- 16) 


The  actual  parameter  vector  determined  by  Wei  and 
Prater  is: 


al 

10 .344  ~ 

a  = 

a2 

= 

3.724 

_a3  _ 

5.616 

( IV- 17 ) 


Using  the  parameters  given  in  (IV-17)  data  was  generated 
from  seven  different  initial  conditions.  This  data  was 
then  perturbated  with  the  use  of  a  random  number  generator 
in  order  that  the  effects  of  experimental  error  could  also 
be  observed. 

It  is  interesting  to  note  that  this  is  a  linear 
system,  but  when  the  parameter  determination  problem  is  formu¬ 
lated  as  a  boundary  value  problem,  the  problem  is  nonlinear 
because  of  the  cross  terms  between  the  state  vector  and  the 
parameter  vector.  This  means  that  the  Jacobians  for  this 
problem  are  not  constant  matrices  as  shown  below: 


■ 
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J  . 
XD 


Hi  +  a2  (1+ke20  L?lk 


el 


a2ke2 


H 


& 


1 k  1 
1  el 


+  a^  (l+ke3Q 


( IV-18 ) 


j 


=  1, 


2 


.  .  .  7 


J 

*3 


Cke2X2 


‘  xl3  &e2(1"x2  "  xl^  PO 


&1  '  kelx2^  C°I1  &e3  (k  - 


x  -  x, 
2  1 


)  -  xj 


(IV-19) 


j  =  1,  2  ...  7 


Data  were  generated  for  the  seven  different  initial 
conditions  given  in  Table  5.  This  data  was  then  perturbated 
with  a  random  relative  error  with  a  normal  distribution. 

These  errors  were  generated  using  the  IBM  360  scientific  sub¬ 
routine  package  (10)  with  a  mean  of  zero  and  standard  devia¬ 
tions  of  0.001,  0.01,  0.1. 


£S  £'  •*£*  X*- 
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TABLE  5 

Initial  Conditions  for  1-Butene  Problem 


Data  Set  No. 

x1  (0) 

x2(0) 

1 

1.0 

0.0 

2 

0.0 

1.0 

3 

0.0 

0.0 

4 

0.24 

0.76 

5 

0.3492 

0.6508 

6 

0.0 

0.4937 

7 

0.413 

0.0 

The  phase  plane  diagram  shown  in  Figure  8  shows  the  rela¬ 
tive  trajectories  of  x^ (t)  and  (t)  .  The  straight  line 
trajectories  represent  eigenvector  directions. 

Using  the  initial  guess  parameter  vector: 


1.0 


a 


(0) 


1.0 


1.0 


( IV-20 ) 


and  all  the  data  sets,  the  solutions  in  Table  6  were  obtained 
using  quasilinearization  directly  without  the  aid  of  data 
perturbation . 


-49- 


i-butene 


-BUTENE 


Figure  8.  Phase  Plane  Diagram  for  1-Butene  Problem  (These  Trajectories  were  Produced  by  the  True  Parameters.) 
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TABLE  6 

Solutions  Using  All  Available  Data 
for  1-Butene  Problem 


Standard 

Deviation 

No .  of 
Iterations 

ai 

a2 

a3 

Sum 

of  Errors 
Squared 

0.1 

7 

10.930 

3.619 

5.540 

0.09481 

0.01 

6 

10.404 

3.713 

5.609 

0.00095 

0.001 

5 

10.350 

3.722 

5.615 

0.00009 

actual  values 

10.344 

3.724 

5.616 

- 

From  the  results  in  Table  6,  it  can  be  seen  that  when 
all  the  data  is  considered  together.,  the  quasilinearization 
procedure  is  relatively  insensitive  to  experimental  error. 

In  order  to  established  the  relative  importance  of 
each  of  the  data  sets,  several  additional  runs  were  made. 

The  results  of  these  numerical  experiments  verified  the  logical 
conclusions  made  earlier  in  this  section. 

Data  sets  five  and  six  were  generated  from  the  initial 
conditions : 


x5(°) 


0.3492 

0.6508 


-6  (0) 


0.0 


( IV- 21 ) 


and 


0.4937 


■ 


■ 


51 


respectively.  These  initial  conditions  are  in  the  two  dif¬ 
ferent  eigenvector  directions.  Attempts  to  find  the  para¬ 
meter  vector  using  only  one  of  these  two  data  sets  proved 
fruitless,  no  matter  how  close  the  initial  guess.  In  other 
words,  a  solution  could  not  be  found. 

However,  by  using  both  of  these  data  sets,  a  solution 
was  easily  found  as  shown  in  Table  7.  The  initial  guess 
vector  was  the  same  as  before. 

TABLE  7 

Solutions  Using  Data  from  the 
Eigenvector  Directions 


Standard 

Deviations 

No .  of 
Iterations 

*i 

a2 

a3 

Sum 

of  Errors 
Squared 

0.1 

9 

12.227 

2.242 

6.146 

0.02946 

0.01 

6 

10.513 

3.588 

5.658 

0.00030 

0.001 

5 

10.361 

3.710 

5.620 

0.000003 

actual  values 

10.344 

3.724 

5.616 

- 

This  confirms  that  data  must  be  taken  in  a  manner  such 
that  all  the  different  eigenvector  directions  are  adequately 


defined . 


' 
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It  would  seem  possible  that  one  experiment,  if 
properly  chosen,  could  be  used  to  find  a  solution.  This 
experiment  would  have  to  be  chosen  so  that  its  data  set 
contained  sufficiently  strong  components  in  each  of  the 
eigenvector  directions.  Data  sets  one  and  two,  with  the 
initial  conditions: 


x-l  (0) 


1.0 

0.0 


x2(°) 


0.0 


1.0 


(IV- 22) 


satisfied  this  requirement.  Using  the  same  initial  guess 
for  the  parameter  vector,  the  results  listed  in  Tables  8 
and  9  were  obtained  for  data  sets  one  and  two  respectively. 


TABLE  8 


Solutions  Using  Data  Set  One 


Standard 

Deviation 

No .  of 
Iterations 

ai 

a 

2 

a3 

D  UIU 

of  Errors 
Squared 

0.1 

oscillating 

without 

^9.6 

^4 . 9 

^3.2 

%0 .013 

0.01 

convergence 

7 

10.167 

3.932 

5.189 

0.001226 

0.001 

6 

10.324 

3.746 

5.569 

0.0000012 

actual  values 

10.344 

3.724 

5.616 

- 
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TABLE  9 

Solutions  for  Data  Set  Two 


Standard 

Deviation 

No .  of 
Iterations 

R  ~ 

a2 

a3 

Sum 

of  Errors 
Squared 

0.1 

No  Co 

n  v  e  r  g 

e  n  c  e 

0  b  t 

a  i  n  e  d 

0 . 01 

7 

10.473 

4.129 

5.575 

0.0001856 

0 . 001 

6 

10.355 

3.758 

5.614 

0.00000186 

actual 

values 

10.344 

3.724 

5.616 

- 

The  remaining  data  sets  were  all  too  close  to  one  or 
the  other  of  the  two  eigenvector  directions  to  get  any  meaning¬ 
ful  results. 

From  the  results  obtained  in  Tables  8  and  9,  it  is 
evident  that  one  experiment,  if  properly  chosen  can  be  used 
to  find  the  parameters.  However,  it  is  also  evident  that 
when  one  such  experiment  is  used,  that  the  results  as  well  as 
the  convergence  may  be  quite  sensitive  to  any  experimental 
error.  In  general,  it  was  observed  that  the  more  experiments 
that  are  included  in  the  analysis,  the  less  sensitive  is  the 
prcoedure  to  experimental  error. 

As  the  proposed  method  depends  on  having  accurate 
initial  conditions  for  each  experiment,  it  would  seem  advis¬ 
able  to  check  the  sensitivity  of  the  results  to  this  assumption. 
To  do  this,  the  first  three  data  sets  were  used  with  exact  and 
perturbated  initial  conditions.  The  perturbated  initial 


. 
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conditions  were 


x  (0) 


1.0  -  y 
0.0  +  y 


-2  <0) 


0.0  +  Y 

1.0  -  Y 


and  (0) 


0.0  +  y 
0.0  +  y 


The  results  for  a  standard  deviation  of  0.001  and  0.1  are 
presented  in  Tables  10  and  11. 


TABLE  10 

Solution  for  1-Butene  Problem  Using  Perturbated 
Initial  Conditions  (a  =  0.001) 


Y 

ai 

a2 

a3 

0.00 

10.346 

3.727 

5.631 

0.01 

10.177 

3.724 

5.545 

0.10 

8.487 

3.968 

4 . 876 

TABLE  11 

Solutions  for  1-Butene  Problem  Using  Perturbated 
Initial  Conditions  (a  =  0.1) 


Y 

ai 

a2 

a3 

0.00 

10.638 

3.930 

5.471 

0.01 

10.479 

3.921 

5.389 

0.10 

8.866 

4.143 

4.749 

. 
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From  these  results  it  is  evident  that  small  inaccur¬ 
acies  can  be  tolerated.  However,  if  a  large  discrepancy 
exists,  the  solution  may  be  strongly  affected. 

In  conclusion,  it  may  be  said,  that  the  design  and 
execution  of  the  experiments  may  determine  whether  or  not  a 
unique  solution  can  be  found.  Further,  if  a  solution  is 
found,  the  dependability  of  the  results  will  be  reflected  by 
the  design  of  the  experiment.  Experimental  design  is  an  ex¬ 
tensive  topic  in  itself  and  it  will  not  be  discussed  further 
here . 


C.  A  Non-Isothermal  Reaction 


The  two  previous  examples  were  somewhat  artificial  in 
that  actual  experimental  data  was  not  used.  The  example 
presented  in  this  section  uses  data  for  the  non-isothermal 
pyrolysis  of  propane  obtained  by  Kershenbaum  and  Martin 
(14,  15).  Thev  suggested  the  following  rate  expression- 


df 

dL 


S 

F 


A ' 


exp 


RT (L ) ' 


P(L) 

a 

1  -  f 

R'T (L) 

1  +  N0  +  f 
F 

P  (L ) 


P  . 
l 


( IV- 23 ) 
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f 

L 


S 


T  (L)  = 

P(L)  = 

P  . 
l 

P 

0 

F 


R 


R' 


A  ' 

E 

a 


fractional  conversion 

distance  from  reactor  entrance  (cm) 

total  length  of  reactor  (cm) 

2 

reactor  cross-sectional  area  (cm  ) 

temperature  along  reactor  (°K) 

pressure  along  reactor  (mmHg) 

inlet  pressure  (mmHg) 

exit  pressure  (mmHg) 

feed  rate  of  propane  (g-moles/sec) 

feed  rate  of  inerts  (g-moles/sec) 

_2  Real 

ideal  gas  constant  (1.987  x  10  -  ) 

(g-mole) (°K) 

ideal  gas  constant  (62361  ) 

(g-mole) (°K) 

.  1-a  -1 

kinetic  rate  constant  (g-mole/cc)  sec 

activation  energy  (Kcal/g-mole) 

order  of  reaction 


The  temperature  was  measured  at  twenty-seven  equally 
spaced  points  along  the  reactor;  T(L)  could  thus  be  accurately 
established.  Conversion  was  measured  only  at  the  outlet. 

Data  for  sixteen  different  experiments  were  used  in  the  analy¬ 
sis  . 

Using  the  notation  given  before,  the  state  vector  is 


x 


f 


(IV- 2 4 ) 


- 

. 

■ 

. 

. 

9  X 
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and  the  parameter  vector  is 


A' 


a 


E 


a 


(IV- 25) 


When  the  problem  is  formulated  as  in  equation 
(IV-23) ,  certain  computational  problems  arise.  These  prob¬ 
lems  were  also  reported  by  Lee  (19)  who  suggested  that 
double  precision  arithmetic  is  needed.  However,  it  was 
found  that  if  A'  and  E  are  normalized  to  the  same  approxi¬ 
mate  order  of  magnitude,  that  higher  precision  was  not 
necessary.  This  was  done  by  letting 


A '  =  10  ^  exp (A) 


( IV- 26 ) 


and  substituting  into  equation  (IV-23) 


a 


df 

dL 


FxlO 


exp (A  - 


E 


RT  (L) 


" P(L) ' 
R'T (L) 


1  -  f 


—i  a 


1  No  +  f 


L 


i 

J 


(IV- 27) 

With  this  formulation,  no  computational  problems  were 
encountered.  The  parameter  vector  now  is 


a  = 


A 

E 


(IV-28) 


a 
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This  problem  was  broken  into  two  parts  by  first  con¬ 
sidering  it  with  the  order  of  the  reaction  fixed  at  unity, 
and  then  considering  the  effects  of  allowing  the  third 
parameter  a  to  be  free. 

With  the  first  problem,  the  unknown  parameter  vector 
is  reduced  to 


( IV-29 ) 


Several  indirect  methods  were  used  by  Kershenbaum 
and  Martin  (14,  15)  in  attempts  to  find  suitable  values  for 
A'  and  E.  All  of  these  methods  took  advantage  of  the  fact 
that  equation  (IV-23)  is  separable  and  an  explicit  expression 
for  A'  in  terms  of  E  could  be  found.  This  yielded  a  set  of 
algebraic  expressions  of  the  form 


A'  =  <|>i  (E) 

i  =  1 ,  2  ...  n 


( IV-30 ) 


with  one  equation  being  for  each  of  the  n  different  experi¬ 
ments  . 

An  attempt  to  find  A'  and  E  by  minimizing 
n 

Z  =  z 
i=l 

by  using  steepest  descent  required  an  excessive  amount  of 
digital  computer  time  and  this  method  was  abandoned. 


[♦±  (E>  -  A'  j 


(IV-31) 


.  - 


■■ 
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They  finally  resorted  to  a  method  which  took  advan¬ 
tage  of  the  approximate  relation  given  below 


log  4>±  (E)  K  viE  +  ui  ( IV-32 ) 

i=l ,  2  .  .  .  n 

A  straight  line  could  be  plotted  for  each  experiment. 
The  slope  and  intercept  could  be  measured  and  then  A'  and  E 
could  be  determined  by  a  linear  least  squares  analysis. 

This  ingenious  approach  will  give  a  reasonable 
approximation  to  A'  and  E.  However,  it  will  not  necessarily 
give  the  best  unbiased  estimate  of  the  parameters  in  a 
least  squares  sense. 

Using  the  algorithm  proposed  in  this  work,  a  direct 
solution  to  this  problem  can  be  obtained.  With  the  least 
squares  criterion,  the  objective  function  to  be  minimized  is 


n 

E 

i=l 


f .  <LJ 

1  m 


-t  2 


f  •  (L  ) 
l  m 


(IV-33) 


with  f ^ (Lm)  and  f^ (Lm)  being  the  predicted  and  experimental 
conversions  for  the  reactor. 

The  solution  obtained  was 


A' 

6.3  x  108 

_  E 

1 

CM 

• 

( IV- 3 4 ) 


with  Z  =  0.034 


i  ■ 


. 
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This  is  somewhat  better  than  the  solution  obtained  by  the 
indirect  approach  which  was 


A ' 

E 

with  Z  =  0.27 

Quasilinearization  combined  with  data  perturbation 
was  used  to  obtain  the  solution  given  in  equation  (IV-34) . 
The  same  solution  was  obtained  with  several  different  ini¬ 
tial  guesses. 


2.4  x  10 


52.1 


11 


( IV- 35 ) 


TABLE  12 

Convergence  Rates  for  the  First 
Propane  Problem 


Initial  Guess 

A 

E/R 

No.  of  Steps 

IBM  360/67 
time  (min) 

35 

26 

8 

3.90 

10 

10 

17 

4 . 56 

30 

30 

16 

4 . 34 

40 

40 

16 

4.44 

0 

0 

36 

7.73 

50 

50 

16 

4.03 

35.40 

26.22 

15 

3.98 
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From  Table  12,  it  is  evident  that  the  solution  can 
be  obtained  in  a  reasonable  amount  of  computational  time. 
However,  it  is  not  a  trivial  problem.  When  the  third  para¬ 
meter  a  was  allowed  to  be  free,  the  computational  task 
becomes  larger. 

The  solution  obtained  for  the  second  problem  was 


A  ' 

“2.13  x  109  ~ 

E 

= 

42.68 

a 

_ 1 . 109 

with  Z  =  0.0303 

This  solution  was  obtained  from  each  of  the  initial  guess 
listed  in  Table  13. 


TABLE  13 


Convergence  Rates  for  the  Second 
Propane  Problem 

Initial  Guess 


8.35 
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It  was  found  that  if  the  initial  guess  for  a  was  greater 
than  1.5,  a  very  small  step  size  was  required  and  thus 
an  excessive  amount  of  computational  time  was  needed. 

Allowing  a  to  be  free  yielded  only  a  small 
improvement  in  the  objective  function.  From  this,  it 
could  be  concluded  that  the  first  order  model  is  suffi¬ 
cient  to  represent  the  data  from  the  given  sixteen  experi¬ 
ments  . 

From  this  example,  it  is  evident  that  reasonably 
small  non-isothermal  kinetic  systems  can  be  handled  with 
the  proposed  algorithm.  For  larger  systems,  with  several 
different  unknown  activation  energies,  a  large  amount  of 
computational  time  may  be  required  due  to  the  extreme  be¬ 
havior  of  the  exponential  function. 

As  the  computational  time  is  directly  proportional 
to  the  number  of  experiments  considered,  it  would  seen  de¬ 
sirable  to  use  the  minimum  number  of  experiments  needed 
to  define  the  parameter  vector  uniquely.  Then  once  an  approxi¬ 
mate  solution  is  obtained,  the  additional  experiments  could  be 

used  to  improve  the  initial  result. 

D .  A  Catalytic  Reaction 

Solid  catalyzed  gas  reactions  such  as: 

-  ci2  +  h2o 


2HC1  +  h  0 


( IV- 3 7 ) 


. 

■ 


63 


in  the  presence  of  a  chromic  oxide  catalyst,  normally 
can  be  modeled  with  a  Langmuir-Hinshelwood  type  rate 
expressions  such  as  the  equation 


P  P 
h20  Cl 


Ke  PHC1 


) 


P  P' 
H,0  CL 

1  +  a2PHCl  +  a3PH  °  +  a4  “  ~2 

2  e  HC1 


( IV-38 ) 


P 


h2o 


P 

HC1 


=  rate  of  reaction  (g-mole  of  oxygen/g-catalyst) 
/ (min) 

=  known  equilibrium  constant  (atm * 1  2) 

=  partial  pressure  of  water  (atm) 

=  partial  pressure  of  oxygen  (atm) 

=  partial  pressure  of  chlorine  (atm) 

=  partial  pressure  of  hydrogen  chloride  (atm) 


Equation  (IV-38)  was  proposed  by  Jones,  Bliss  and  Walker 
(11)  to  model  the  above  reaction  with  a  fourth  order  un¬ 
known  parameter  vector. 


a 


(IV-39) 


. 
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For  this  type  of  problem,  although  the  nonlinear 
differential  equation  modeling  the  system  is  of  first 
order  the  parameter  vector  is  of  high  order.  Models  with 
up  to  seven  constants  are  not  unusual.  An  extensive 
experimental  and  analytical  research  program  is  required 
if  such  a  model  is  to  be  identified  properly.  For  this 
reason,  Levenspiel  (21)  suggests  that  simpler  empirical 
expressions  could  be  satisfactorily  fitted  to  the  data 
for  the  purpose  of  design.  However,  if  the  data  is  to 
be  extrapolated,  a  mechanistic  model  is  desirable. 

Extrapolation  is  at  best  a  dangerous  procedure. 
But,  a  mechanistic  equation  does  give  some  qualitative 
insight  into  what  may  be  expected  under  other  conditions, 
and  may  be  used  to  design  further  experimental  studies. 

For  these  reasons,  the  application  of  the  proposed  algo¬ 
rithm  to  this  type  of  eauation  was  considered. 

The  analysis  of  this  type  of  problem  for  the  most 
part  has  been  limited  to  the  least  squares  analysis  of 
algebraic  relations  similar  to  equation  (IV-38) .  The 
rates  are  measured  and  the  parameters  are  determed  by  non¬ 
linear  least  squares.  Frequently,  the  relations  may  be 
manipulated  to  form  a  linear  system  with  respect  to  the 
parameters.  Rewriting  equation  (IV-38)  as 


1 

> 

- 

J 
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l/a1  (1  +  a2PHCl  +  a3PH20  +  a4 


P  P 
H2°  C12 

K  P^ 
e  HC1 


) 


P  P 

1  H.O^Cl 

- -  (P2  -  - - - - 

v  0  2 
2  K  P , 


(IV-40) 


0. 


e  HC1 


produces  such  a  linear  system,  but  with  the  sacrifice  of 
the  desired  fitting  criterion. 

If  the  rates  cannot  be  measured  directly-  they 
must  be  obtained  by  numerically  differentiating  the  inte¬ 
grated  data,  thus  introducing  a  large  amount  of  error. 

Under  these  circumstances ,  it  would  seem  desirable  to 
attack  the  integrated  data  directly.  Such  an  approach 
leads  to  the  inevitable  boundary  value  problem  to  which 
CTuasilinearization  and  data  perturbation  may  be  applied. 

As  an  example,  the  problem  presented  in  the  ecxuation 
(IV-37)  and  (IV-39)  was  considered.  Taking  the  partial 
pressure  of  oxygen  to  be  the  dependent  variable 


with 


NB  d  p02 

W'P.  dt 

1 


( IV-41 ) 


N  =  total  number  of  gm-moles  initially  present 
B 

in  the  reactor 

W'  =  weight  of  catalyst  present  in  the  reactor  (gm) 

P.  =  the  initial  pressure  in  the  reactor  (atm) 
i 


' 
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If  the  system  were  free  from  external  interference,, 

the  total  pressure  and  the  partial  pressure  of  each  of  the 

other  components  in  the  reactor  would  be  a  function  of 

Pn  *  Presuminq  this,  the  total  pressure  would  be  given 
u2 
by 


P  . 
1 


tot 


N 


B 


N 


”b  -  (No,  -  P1  po  > 

-  2  i  2  — 1 


(IV-42) 


However,  the  pressure  measurements  throughout  the  experi¬ 
ment  did  not  follow  this  relation  too  closely,  indicating 
some  external  inter ference,  probably  stemming  from  the 
operation  of  the  recirculation  pump.  It  was  found  that 
the  pressure  data  could  be  better  represented  by  a  rela¬ 
tion  of  the  form 


P 


tot 


0 


2  -J 


B  ( t ) 


( IV- 4  3 ) 


where  B(t)  is  a  correction  factor  correlated  from  the 
actual  pressure  data.  The  correction  factor  never 
exceeded  15%  and  was  correlated  to  be  within  one  and  a 
half  percent  with  polynomials  of  up  to  the  third  order. 
The  error  in  the  correction  factors  was  small  compared 
with  the  experimental  error  in  the  data. 

The  partial  pressures  of  the  other  components  are 


given  by  the  relations 
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P 

HCl 


P 


Cl 


2 


P 


h2o 


N 

HCl 

,  >»0j 

nb 

P  . 

1 

V 

I 

JQ 

f— • 

(U°2  - 

P  . 

1 

V 

Nci  + 

2 

2  1"»2 

nb 

P  . 

1 

V 

hb  - 

<H°2  - 

nb 

p . 

1 

V 

Nh2o  + 

2  (N°2 

nb 

p . 

1 

V 

nb  - 

‘S  - 

nb 

p . 

1 

V 

P.  .  ( IV-44 ) 

tot 


Ptot  (IV-45) 


Ptot  (IV'46) 


with  N 
N 


°2 

HCl 


U 


H2° 


no.  of  moles  of  02  initially  in  the  reactor 

no.  of  moles  of  HCl  initiallv  in  the  reactor 

no.  of  moles  of  CL2  initially  in  the  reactor 

no.  of  moles  of  H20  initially  in  the  reactor 


The 

reactor  is 


total  number  of  moles  initially  present  in  the 
given  by 


nb  N02  +  lJHC1  +  NC1  +  IJH20  +  NI 


where  N  =  the  number  of  moles  of  inerts  initially  present 


in  the  reactor. 
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Combining  eauations  (IV-44),  (IV-45),  and 

with  (IV-43)  produces 


HCl 


nhci  4 (wo. 


p . 
1 


v 


%Biti 


Cl. 


N 


Cl. 


2  (N 


0. 


p . 
1 


p0  } 
2 


P  . 


N 


-  B  (t) 


B 


h2o 


1JH20  +  2  (N02 


p . 
1 


p0  > 

2 


N 


B  ( t ) 


B 


( IV-46 ) 


( IV- 4 7 ) 


(IV- 48) 


(IV- 49) 


Partial  pressure  data  was  made  available  for  nine  experi¬ 
ments  at  355°C,  eleven  experiments  at  340°C,  and  nine 
experiments  at  325°C.* 

Jones,  Bliss  and  Walker  (11)  determined  the  rate 
parameters  for  each  temperature  by  using  linear  least 
squares  with  equation  (IV-40).  The  rates  were  obtained  by 
plotting  the  data  and  estimating  the  slope.  Subsequently 
Mezaki  and  Bliss  (23)  fitted  the  data  for  all  three  tem¬ 
peratures  using  nonlinear  least  sauares. 

Both  of  these  correlations  made  use  of  numerical 
differentiation  of  the  observed  data  and  therefore  suffer¬ 
ed  from  the  inevitable  loss  of  accuracy. 


* 


Data  was  taken  by  A  M.  Jones  and  was  kindly  supplied 
by  Professor  Harding  Bliss  of  Yale  University. 


< 
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Attempts  to  solve  for  the  parameters  for  each 
temperature  separately  using  quasilinearization  and 
data  perturbation  yielded  mixed  results . 

For  the  data  at  355°C,  the  following  results 
were  obtained. 


a 


9.8  x  10 
13 . 2 


47.9 

762 


(IV- 50) 


with  a  least  scruare  error  of  0.000747  .  The  parameter 
vector  reported  by  Jones,  Bliss,  and  Walker  (11)  was 


a 


2.4  x  10 

3.1 

25 

230 


( IV-51 ) 


with  a  least  scruare  error  of  0.00315. 

Using  the  data  given  for  340°C,  the  results  were 


a  = 


5.57  x  10 
0.791 
121 . 3 


-4 


( IV-52 ) 


-34.8 


. 
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with  a  least  square  error  of  0.00144.  Jones,  Bliss  and 
Walker  (11)  found 


a  = 


0.92  x  10 

0.22 


21 

142 


(IV- 53) 


with  a  least  scruare  error  of  0.00708.  The  results  given 
in  ecruation  (IV-52)  contain  a  negative  rate  constant. 

This  is  inconsistent  with  kinetic  theory.  However,  given 
this  set  of  data,  this  model  and  the  unweighted  least 
square  criteria,  the  parameter  vector  given  in  ecruation 
(IV-52)  is  an  optimum  for  reproducing  the  experimental 
data . 

For  the  data  given  at  325°C,  the  parameters  tended 
to  grow,  seemingly  without  bound,  with  each  successive 
step.  Using  the  parameters  found  by  Jones,  Bliss,  and 
Walker  (11)  as  an  initial  guess,  the  results  of  several 
different  steps  are  presented  in  Table  14. 
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TABLE  14 

Hydrogen  Chloride  Problem  Results  for  325°C  Data 


Step 

No. 

ai 

a2 

a3 

a4 

Sum  of 
Errors 
Squared 

-4 

0 

0.26x10 

0.01 

11.0 

18  o  0 

0.0183 

4 

0 .68xl0"4 

1.71 

24 . 5 

672 

0.0135 

5 

3 . llxlO-4 

10.7 

101 

2,259 

0 . 0097 

-4 

7 

34x10 

124 

1076 

19,964 

0.0083 

It  can  be  seen  from  Table  14,  that  there  are  at 
least  several  parameter  vectors  which  will  reproduce  the 
data  better  by  the  least  square  criteria  than  those  given 
by  Jones,  Bliss,  and  Walker  (11).  However,  most  of  these 
other  parameter  vectors  are  physically  inconsistent  with 
the  results  obtained  at  the  other  temperatures. 

From  this,  it  appears  that  the  integrated  rate 
data  is  taken  over  a  domain  which  is  insufficient  to  excite 
all  modes  of  the  system.  A  similar  situation  may  exist 
with  the  data  given  at  340°C,  producing  the  inconsistent 
negative  rate  constant.  If  this  is  so,  the  parameters  can¬ 
not  be  properly  identified. 

The  convergence  for  the  data  at  355°C  and  340°C  was 
quite  rapid  when  an  order  of  magnitude  guess  was  supplied. 
There  is  then  good  reason  to  suppose  that  the  proposed 
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procedure  will  be  quite  useful  in  the  analysis  of 
Langmuir-Hinshelwood  models  when  integrated  rate  data 
is  obtained.  However,  the  discouraging  results  obtained 
from  the  data  given  at  325°C  and  340°C  indicates  that  a 
considerable  amount  of  caution  should  be  exercised  when 
attempting  to  correlate  with  this  type  of  model. 

In  principle,  it  would  not  be  too  difficult  to 
modify  the  proposed  procedure  so  that  constraints  on  the 
parameter  vector  could  be  incorporated.  If  this  were 
done,  rates  obtained  by  direct  measurement  could  be 
incorporated  so  that  all  integrated  and  differential  data 
could  be  fitted  by  the  same  correlation  technique.  This 
would  be  of  considerable  advantage  in  attacking  problems 
such  as  the  one  above . 

E .  A  Complex  Nonlinear  Reaction 

To  identify  the  parameters  in  a  complex  reaction 
model  is  at  best  a  very  difficult  problem.  This  will  be 
illustrated  by  the  problems  encountered  in  attempts  to 
identify  a  system  with  the  proposed  mechanistic  model. 


. 

' 
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H+  +  x. 


a, 


x 


W  +  x 


X2  +  x5 


( IV- 5 4 ) 


x, 


x6  +  X2  +  H 


+ 


x, 


x  -  +  W  +  H 


+ 


Several  experiments  were  performed  at  each  of  the 
two  temperatures  80°C  and  100°C  in  a  batch  reactor.  The 
concentrations  of  x^ ,  *2 ,  and  x^  were  measured  at  several 
different  times  during  the  reaction.  The  components  x^ 
and  Xj-  were  unmeasurable  intermediates  and  x^  was  given 
by  the  following  linear  relation  derived  from  the  stoichio¬ 
metry 


x,  =  C  (x,  +  x-  +  x.  +  x,-)  (IV-5  5) 

6  o  1  3  4  d 

Hydrogen  ion  concentration  (H+) ,  water  concentration  (W) 
and  the  constant  Cq  remained  fixed  throughout  any  one 
experiment.  At  the  request  of  Chemcell,  who  generously 
supplied  the  data  for  this  example,  the  other  components 
will  not  be  identified. 
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From  the  proposed  mechanistics  model,  the  system 
of  ordinary  differential  equations  may  be  written* 


dx-^ 

dt 


dx2 

dt 


dx3 

dt 


dx 


4 


dt 


dx 


5 


dt 


-  a.,H  xx 


a5  X4 


a2Wx4 


a3X5 


"  a6X2X5 


a7H+X6X2 


x 


5 


alH+Xl 


arx,  -  a„Wx,  +  a,x„xr 
5  4  2  4  6  2  5 


a2Wx4  ~  a6X3X5  ~  a3X5  +  a7H  X6X2  ~  a4x5 


(IV-56) 


(IV-57) 


( IV- 5 8 ) 


( IV-59 ) 


( IV- 6  0 ) 


with  the  unknown  parameter  vector  being 


a 


(IV-61) 
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and  the  state  vector  being: 


x  = 


( IV- 6 2 ) 


Attempts  to  solve  this  system  directly  met  with 
difficulties,  as  the  system  of  ordinary  differential 
equations  was  unstable  to  numerical  integration.  This 
was  caused  by  the  tremendous  difference  in  magnitude  of 
the  time  constants  associated  with  the  different  com¬ 
pounds.  To  avoid  this  problem,  the  stationary  state 
approximation  was  employed  by  setting 


dx 


4 


dt 


0 


( IV-63 ) 


and 


0 


( IV- 6 4 ) 


and  using  the  resultant  system  of  algebraic  equations  to 
obtain  the  following  relations. 
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(a7H+x2  -  a2W)a^H+x^  -  a7H+x2 (CQ  -  x^  -  x^) (a^  +  a2) 
(a2W  “  a7H+x2 ) agx2  ~ (agx2  +  a3  +  a4 +  a7H+x2 )  (a^  +  a2W) 


al^  X1  a6X2 

ac  +  a~W  aP  +  a«W 

o  Z  5  2 


( IV- 6  5 ) 

(IV- 6 6 ) 


Substituting  (IV-65)  and  (IV-66)  into  the  differen¬ 
tial  equations  (IV-56) ,  (IV-57),  and  (IV-58)  produced  a 

stable  system  with  the  state  vector  being  reduced  to 
order  three 


(IV- 6 7 ) 


Although,  with  this  modification,  no  integration 
stability  problems  were  encountered,  attempts  to  find 
all  seven  parameters  were  not  successful.  If  the  data 
sets  are  examined,  it  will  be  found  that  all  the  initial 
conditions  are  of  the  form- 


x(0)  = 


( IV- 6 8 ) 


0 


. 
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It  was  suspected  that  data  taken  from  this  initial  con¬ 
dition  alone  was  not  sufficient  to  identify  distinctly 
the  various  forward  and  backward  rate  parameters  in  the 
first  two  steps  in  the  mechanistic  model. 

A  simplified  model  was  then  proposed* 


+  1 

x1  +  H  +  W  — =->  x2  +  x4 


-4  *2  +  x6  +  H' 


3  + 

x4  ^  x3  +  W  +  H 


which  may  be  described  mathematically  by  the  differential 
equations : 


dx. 


dt 


-  axH  Wxx 


( IV- 6  9 ) 


dx. 


dt 


a2X4  ”  a5H+x2X6  +  a^H+Wx]_ 


( IV- 7  0 ) 


dx. 


dt 


a,x 
3  4 


( IV- 71 ) 


and  the  algebraic  relations 


X£  =  c~ 
6  o 


(xl  +  x3 


+  x4) 


( IV- 7 2 ) 
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(IV-73) 


The  relation  (IV-72)  comes  from  the  stoichiometry  and  the 
expression  for  x^  (IV-73)  again  comes  from  the  stationary 
state  approximation.  The  parameter  vector  has  now  been 
reduced  to 


a 


1 


a 


2 


( IV-74 ) 


a 


a 


3 


a 


5 


In  attempts  to  identify  the  parameters  in  the  vec¬ 


tor  (IV-74),  an  interesting  phenomenon  was  encountered. 

The  parameters  a^  and  a,.  were  accurately  identified  for 
the  data  given  at  the  two  different  temperatures.  For 
the  parameters  a-2  and  a^  a  large  number  of  different  solu¬ 
tions  could  be  obtained.  Each  yielded  the  same  least 
squares  minimum. 

On  closer  examination  of  the  model,  it  will  be 
found  that  the  following  relation  is  approximately  true 

1 


x 


4 


( IV- 7 5 ) 


a 
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and  only  appears  in  two  terms  where  is  is  dominant. 
The  first  is  the  term  i n  equation  (IV-70)  and  the 
second  is  a^x^  in  equation  (IV-71) . 

From  the  above; 


a2X4 


and 


a3x4 


a 


a2  +  a3 


a 


a2  +  a3 


( IV-76 ) 


( IV-77 ) 


Since  the  concentration  of  x^  is  not  specified, 
other  than  being  several  orders  of  magnitude  smaller  than 
the  principal  components,  as  long  as  the  terms  on  the 
right  hand  sides  of  equations  (IV-76)  and  (IV-77)  remain 
constant,  the  model  will  reproduce  the  data  equally  well. 
The  constant  a^  was  found  to  be  much  larger  than  83. 
Therefore  the  approximations 


a 


2 


a 


2 


+ 


1 


(IV-78) 


and 


(IV- 79) 
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may  be  made.  So,  if  the  ratio  of  a^  to  is  the  samgr 
different  parameter  vectors  will  produce  the  same  solution 
profiles.  This  is  what  was  observed. 

To  verify  these  observations,  examine  the' 
following  table  of  solutions  for  the  data  at  100°C. 

The  concentration  profiles  produced  by  any  one  of  the 
three  parameter  vectors  given  in  Table  15  are  shown  in 
Figures  9,  10,  11,  12,  13  and  14. 

TABLE  15 

Chemcell  Problem 

Various  Solutions  for  100°C  Data 

Sum  of 


al 

a2 

a3 

a5 

a3/a2 

Errors 

Sauared 

0.0815 

0.982 

0.151 

7.15 

0.154 

0.00023 

0.0815 

2.11 

0.323 

7.11 

0.153 

0.00023 

0.0815 

0.174 

0.027 

7.47 

0.155 

0.00023 

Similarly,  the  results  for  the  data  given  at  80°C 


verified  the  above  as  shown  in  Table  16. 
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Figure  9  Solution  to  Chemcell  problem  100°C  (data  set  one 
[H+]  =  0.012  g-mole/l  CG=  0.0759  g-mole/l 
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[H+]  -  0.00241  g-mole/l  CQ=  0.0774  g-mole/l 
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Figure  11  Solution  to  Chemcell  problem  100°C  (data  set  three) 
[H+]  =  0.00614  g-mole/l  CG=  0.0828  g-mole/l' 
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Figure  12  Solution  to  Chemcell  problem  100°C  (data  set  four) 
[H+]  -- 0.00986  g-mole/l  CG  =  0.0784  g-mole/l 
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Figurel3  Solution  to  Chemcell  problem  100°C  (data  set  five) 
[H+]  =  0.0125  g-mole/l  CQ=  0.0749  g. mole/I 
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Figure  14  Solufion  to  Chemcell  problem  100°C  (data  set  six 
[H+]  =  0.0103  g-mole/l  Co=  0.01563  g-mole/l 
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TABLE  16 

Chemcell  Problem 
Various  Solutions  for  80°C  Data 


Sum  of 


al 

a2 

a3 

a5 

a3/a2 

Errors 

Squared 

0.0106 

0.183 

0.025 

1.37 

0.136 

0.000072 

0.0106 

1.02 

0.138 

1.30 

0.135 

0.000070 

0.0106 

2.23 

0.302 

1.29 

0.135 

0.000070 

All  of  the  parameters  vectors  given  in  Tables 
15  and  16  reproduced  the  data  equally  well.  When  the 
profiles  were  plotted  for  each  experiment,  the  different 
curves  could  not  be  distinguished. 

For  each  experiment  there  were  data  points  taken 
at  24  and  48  hours.  These  points  were  not  used  in  the 
fitting  procedure  in  order  to  eliminate  the  large  amount 
of  integration  time  required.  However,  once  the  para¬ 
meters  have  been  determined,  thes.e  points  may  be  used  to 
check  the  results. 

With  each  of  the  sets  of  parameters  given  in 
Table  14  for  the  80°C  data,  the  predicted  values  of  the 
data  at  24  and  48  hours  were  the  same.  This  is  shown  in 
Figures  15,  16,  17,  18  and  19. 
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Figure  15  Solution  to  Chemcell  problem  80°C  (data  set  one 
[H+]  =  0.0127  g- mole/I  CQ  =  0.0791  g-mole/l 
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Figure  16  Solution  to  Chemcell  problem  80°C  (data  set  two) 
[H+]  =  0.0230  g-mole/l  CQ=  0.0795  g-mole/l 
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Figure  17  Solution  to  Chemcell  problem  80°C  (data  set  three) 
[H+]  =  0.0368  g-  mole  / 1  CQ  =  0.0822  g-mole/l 
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Figure  18  Solution  to  Chemcell  problem  80°C  (data  set  four) 
[H+]  :  0.0345  g-mole/l  CQ=  0.0825  g-mole/l 
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Figurel9  Solution  to  Chemcell  problem  80°C  (data  set  five) 
[H+]  =  0.0345  g-mole/l  CQ  =  0.0825  g-mole/l 
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Even  though  and  could  not  be  uniquely 
determined,  the  model  has  been  successfully  fitted  to 
the  given  data,  and  may  be  used  to  produce  information, 
with  reasonable  certainty,  in  the  range  of  conditions 
under  which  the  experiments  were  performed. 

To  establish  a^  and  a^  uniquely  measurements 
need  to  be  made  of  the  intermediate  x^  .  This,  however, 
is  not  necessary  if  the  model  is  only  to  be  used  for  the 
above  purposes. 

The  most  important  conclusion  that  can  be  drawn 
from  this  example,  is  that  identifying  complex  systems 
is  not  a  routine  operation.  The  analysis  must  be  closely 
coupled  with  the  design  and  execution  of  the  experiments. 
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V.  IDENTIFICATION  OF  SIMPLE  CONTROL  MODELS  FOR 

COMPLEX  SYSTEMS 


For  many  complex  systems,  the  sophisticated 
mathematical  models  that  are  needed  to  accurately 
simulate  them  recruire  a  large  amount  of  time  for 
solution,  or  they  may  defy  solution  altogether.  Time 
is  not  of  critical  importance  in  most  design  problems 
as  long  as  a  solution  is  obtained,  and  the  cost  of  ob¬ 
taining  that  solution  can  be  justified  economically. 
However,  for  control  problems  with  the  independent 
variable  being  time  itself,  it  is  essential  that  solu¬ 
tions  be  obtained  quickly  so  that  the  appropriate  mani¬ 
pulations  of  the  control  variables  can  be  made  at  the 
proper  times.  For  these  complex  systems  then,  it  is 
important  to  develop  simple  theoretical  or  correlative 
models  which  can  be  used  over  a  narrow  range  of  operating 
conditions . 

In  many  cases,  the  initial  conditions  are  not 
known  for  some  of  the  state  variables.  This  being  the 
case,  the  proposed  algorithm,  in  its  simplified  form, 
cannot  be  used.  However,  there  are  some  cases  for  which 
this  is  not  true.  An  example  is  given  and  discussed 


in  this  section. 
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A .  Flow  of  a  Tracer  Through  a  Packed  Bed 


The  analysis  of  flow  through  packed  beds  is  of 
considerable  importance  in  a  number  of  systems,  the 
least  of  which  is  certainly  not  the  non-homogeneous 
catalytic  reactor. 

The  partial  differential  equation 

2 

9x  1  9  x  9x 

— —  -  -  2  -  (V-l) 

9 1  P  9L  9L 

e 


dimensionless  concentration 
dimensionless  length 
dimensionless  time 
Peclet  number 


e 

is  a  mathematical  model  which  can  be  used  to  describe 
the  unsteady  state  behavior  of  a  tracer  flowing  through 
a  packed  bed.  This  model  has  already  been  simplified. 
The  radial  diffusion  terms  have  been  dropped,  leaving 
only  the  axial  diffusion  and  convective  terms  on  the 
right  hand  side  of  the  equation. 

Equation  (V-l)  can  be  solved  numerically 
without  a  great  deal  of  difficulty.  With  packed  bed 
reactors,  however,  the  nonlinear  rate  expressions 
must  be  added,  and  the  resultant  nonlinear  partial 
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differential  ecruations  are  quite  difficult  to  handle. 

This  example  was  considered  with  the  application  to 
more  complex  svstems  being  the  ultimate  goal. 

Data  was  generated  by  H.S.  Koonar  (17)  for 
the  transient  behavior  of  (V-l)  with  the  initial  con¬ 
dition 

x (0, L)  =  0.75  ( V—  2 ) 

and  the  boundary  conditions 

3x 

—  (t,l)  =  0  (V-3 ) 

9  L 

x  ( t ,0 )  =  f(t) 

=  0.75  +  0.25  SIN  (cot)  (V-4) 

It  was  solved  by  using  central  finite  difference  rela¬ 
tions  to  approximate  the  terms  on  the  right  hand  side 
of  the  equation.  Doing  this  results  in  a  system  of 
linear  ordinary  differential  equations  with  time  being 
the  independent  variable.  Eleven  grid  points  were  used 
and  the  system  of  linear  ordinary  differential  equations 
was  solved  by  finding  the  eigenvalues  and  eigenvectors  of 
the  coefficient  matrix.  In  Table  17  are  the  outlet  con¬ 
centrations  generated  from  w  =  0.1,  1.0,  5.0,  10.0. 


oney 
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TABLE  17 


Data  for  Tracer  Problem  x(t,l) 


t/ll) 

0.1 

o 

• 

i— 1 

5.0 

10.0 

o 

• 

o 

0.7500 

0.7500 

0.7500 

0.7500 

0.5 

0.7954 

0.7602 

0.7523 

0.7888 

1.0 

0.7523 

0.8250 

0.7599 

0.8423 

1.5 

0.7384 

0.9095 

0.7708 

0.6853 

2.0 

0.7575 

0.9692 

0.7828 

0.8100 

2.5 

0.7737 

0.9781 

0.7951 

0.7404 

3.0 

0.7611 

0.9341 

0.8073 

0.7188 

3.5 

0.7370 

0.8459 

0.8194 

0.8204 

4.0 

0.7354 

0.7351 

0.8314 

0.6784 

4.5 

0.7588 

0.6286 

0.8432 

0.8041 

5.0 

0.7735 

0.5525 

0.8548 

0.7447 
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It  has  been  suggested  that  (26),  (24)  many  com¬ 

plex  systems  can  be  modeled  sufficiently  for  control 
purposes  by  a  second  order  differential  ecruation  with 
a  time  delay  such  as 


d  x (t, 1) 


dt 


2  +  a]_ 


dx (t, 1) 


dt 


+  a2  |x(t,l)  -  f(t-tD)]  =  0 

(V-5 ) 


with  tD  being  the  time  delay.  By  defining 


x1  =  x (t, 1) 


(V-6 ) 


and 


dx (t , 1) 
dt 


(V-7 ) 


equation  (V-5)  can  be  written  as  a  system  of  first  order 
differential  equations 


dx^ 

dt 


x 


2 


dx2 

dt 


alX2  ‘  a2  LX1  '  f(t'tD)]  (V’8) 


with  the  unknown  parameter  vector 


(V-9 ) 
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The  time  delay  tD  cannot  be  treated  directly  by  quasi¬ 
linearization  and  data  perturbation.  For  small  prob¬ 
lems  such  as  this,  the  time  delay  can  be  determined  by 
a  search  procedure.  For  more  complex  systems  with 
many  different  unknown  time  delays,  other  methods  must 
be  employed. 

With  several  different  values  for  tD  the  parameter 
vector  was  determined  with  the  initial  guess  being 


a 


1.0 

1.0 


Table  18  contains  these  values. 


(V-10 ) 


100 


TABLE  18 

Solutions  of  Tracer  Problem 


ai 

a2 

Sum  of  the 
Errors  Squared 

0.100 

4.39 

6.69 

0.003203 

0.150 

5.40 

8.59 

0.000891 

0.170 

5.96 

9.59 

0.000466 

0.180 

6.28 

10.18 

0.000371 

0.187 

6.54 

10.62 

0 . 000350* 

0.190 

6 .65 

10.83 

0.000353 

0.200 

7 . 07 

11.57 

0.000418 

0.210 

7.56 

12.42 

0.000560 

0.230 

10.68 

17.79 

0.001863 

*  Value  for  optimum  time  delay 

The  parameters  obtained  for  each  of  the  given 
estimates  of  the  time  delay  all  reproduce  the  data 
well.  This  can  be  seen  by  observing  Table  19  which 
contains  the  deviations  from  the  actual  solution, 
presented  in  Table  17 ,  for  the  optimum  time  delay 
(tD  =  0.187) . 
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TABLE  19 


Deviation  of  the  Model  from  the  Actual  Values 

for  Tracer  Problem 

|x (t, 1)  -  x (t,l)  | 


t/  0) 

0.1 

1.0 

5.0 

10.0 

0.0 

0.0000 

0.0000 

0.0000 

0.0000 

0.5 

0 . 0007 

0.0019 

0.0015 

0.0026 

1.0 

0.0022 

0.0036 

0.0024 

0.0010 

1.5 

0.0062 

0 . 0035 

0 .0024 

0.0013 

2.0 

0.0042 

0.0036 

0.0028 

0.0048 

2.5 

0.0001 

0.0019 

0.0029 

0.0013 

3.0 

0.0002 

0.0017 

0.0028 

0.0035 

3.5 

0.0041 

0.0017 

0.0028 

0 .0031 

4.0 

0.0074 

0.0022 

0.0029 

0.0013 

4.5 

0.0035 

0.0027 

0.0029 

0.0045 

5.0 

0.0004 

0.0032 

0.0030 

0 . 0013 

These  results  confirm  that  a  second  order  dif¬ 
ferential  equation  with  a  time  delay  can  be  used  as  an 
adeauate  control  model  for  the  flow  of  a  tracer  through 
a  packed  bed . 


. 


■ 
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The  computational  time  required  to  solve  this 
problem  was  quite  small  (less  than  one  minute  of 
360/67  computational  time).  Therefore,  with  small 
problems,  it  is  quite  practical  to  find  the  unknown 
time  delay  by  searching  procedure. 


. 
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VI .  DISCUSSION  AND  CONCLUSIONS 

Quasilinearization  and  data  perturbation  has  been 
shown  to  be  a  powerful  tool  in  the  solution  of  certain 
identification  problems.  Certain  limitations,  however, 
were  evident  from  the  experience  gained  in  analyzing 
the  examples  given  in  this  work. 

A .  Numerical  Problems 

There  are  two  major  numerical  operations  which 
must  be  performed.  Both  of  these  operations  can  be  a 
source  of  difficulties. 

The  first  is  the  integration  of  the  large  number 
of  initial  value  problems  associated  with  each  itera¬ 
tion  of  the  quasilinearization  procedure.  In  some  in¬ 
stances,  these  equations  may  be  unstable  to  numerical 
integration.  If  this  is  the  case,  there  are  three 
alternatives.  The  first  is  to  use  a  small  step  size  which 
will  result  in  consuming  a  large  amount  of  computer  time. 
The  second  alternative  is  to  find  a  better  numerical 
integration  procedure.  Perhaps  one  of  the  multistep 
methods  such  as  Hammings  prediction  corrector  method  (7) 
could  be  used  to  replace  the  fourth  order  Runge  Kutta 
method  (8)  presently  being  used.  The  final  alternative 
is  to  make  use  of  a  hybrid  computer. 


lie 
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The  second  numerical  operation  is  the  solution 
of  the  least  squares  problem  at  the  completion  of  each 
quasilinearization  iteration.  It  is  well-known  that 
this  problem  is  subject  to  ill-conditioning  if  the 
number  of  the  unknown  parameters  and  initial  conditions 
becomes  large  (greater  than  seven) .  This  problem, as 
has  been  mentioned  before,  can  be  avoided  by  using  the 
Chebyshev  criteria  for  fitting  the  data. 

B .  Experimental  Design 

Advances  have  been  made  recently  in  the  theory 
of  the  design  of  experiments.  Especially  noted  is  the 
work  done  by  Box,  Hunter,  and  their  co-workers  (3) ,  (9) , 

and  (16) .  It  is  evident  from  their  work  that  if  a 
mechanistic  or  correlative  model  is  the  ultimate  aim  of 
an  experimental  program,  it  is  essential  to  determine 
the  nature  of  the  model  early  so  that  subsequent  experi¬ 
mentation  can  be  designed  appropriately.  With  the  proper 
design  of  experiments,  at  least  some  of  the  identifica¬ 
tion  failures  will  be  eliminated.  Unfortunately,  most 
of  the  work  that  has  been  done  has  been  associated  with 
the  identification  of  systems  which  can  be  represented  by 
algebraic  relations  or  can  be  reduced  to  this  form. 


' 
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Until  the  theory  of  experimental  design  has 
been  extended  to  function  space,  intuition  and  judge¬ 
ment  must  be  used  in  the  identification  of  systems 
which  can  be  described  adequately  only  by  a  system  of 
ordinary  differential  equations.  The  theory  for  linear 
systems  is  more  advanced  (20),  (27)  and  can  be  used  as 

a  guide  to  the  solution  of  problems  which  are  of  a  non¬ 
linear  nature. 

C .  Computer  Time  Requirements 

Few  of  the  problems  presented  in  this  work  requir¬ 
ed  more  than  five  minutes  of  time  on  the  IBM  360/67  com¬ 
puter.  However,  there  were  indications  that  the  time 
requirements  may  be  disproportionately  larger  for  more 
complex  systems.  If  this  is  true,  then  ways  must  be  devised 
to  improve  the  efficiency  of  this  procedure  or  another  more 
effective  method  must  be  found. 

Two  obvious  ways  of  cutting  down  the  time  require¬ 
ments  are  to  improve  the  programming  and  to  use  more 
efficient  numerical  techniques.  Some  improvement  un¬ 
doubtedly  could  be  made  in  this  direction. 

As  the  amount  of  time  required  is  directly  pro¬ 
portional  to  the  number  of  experiments  being  considered, 
it  would  seem  advisable  to  minimize  the  total  number 


. 
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used  to  identify  the  system.  Once  a  few  experiments 
had  been  used  to  get  an  approximate  result,  all  of 
the  experiments  could  then  be  used  to  find  the  true 
solution.  This  approach  may  be  profitable  in  complex 
problems,  if  the  few  experiments  chosen  are  sufficient 
to  identify  the  system.  The  choice  of  these  few  experi¬ 
ments  is  closely  coupled  to  the  problem  of  experimental 
design  discussed  above. 

The  method  for  adjusting  the  step  size  for  the 
data  perturbation  is  completely  arbitrary.  If  a  technique 
were  developed  for  choosing  the  optimum  step  size  for  each 
step,  a  considerable  advantage  would  be  gained. 

D .  Future  Work 

To  extend  the  present  work  to  more  general  systems 
would  be  of  obvious  advantage. 

In  principle,  the  modification  of  the  method  to 
handle  more  general  boundary  conditions,  both  linear  and 
nonlinear,  is  not  difficult.  However,  the  addition  of 
nonlinear  boundary  conditions  and  constraints  will  pro¬ 
bably  require  larger  amounts  of  computer  time. 

It  has  been  noted  that  the  data  perturbation 
procedure  will  work  satisfactorily  only  if  reasonably 
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large  steps  may  be  taken.  If  the  number  of  steps 
becomes  large,  the  computer  time  required  will  become 
excessive . 

To  circumvent  this  problem,  it  may  prove  advan¬ 
tageous  to  combine  this  method  with  one  of  the  hill 
climbing  methods  such  as  Rosenbrock's  (25).  in  general, 
the  rate  of  convergence  of  the  hill  climbing  procedures 
tends  to  be  slower  as  the  solution  is  approached.  This 
is  quite  the  contrary  to  the  rate  of  convergence  of  the 
quasi linearization  procedure,  which  tends  to  accelerate 
as  the  solution  is  approached.  A  combination  of  these 
two  methods  may  yield  fruitful  results. 
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NOMENCLATURE 

A 

normalized  Arrhenius  rate  constant 

A' 

Arrhenius  pre-exponential  rate  constant 

a 

parameter  vector 

B 

pressure  correction  factor 

C 

max 

maximum  allowable  relative  change  in  the 
parameters  for  one  iteration 

C  . 
mm 

minimum  allowable  relative  change  in  the 
parameters  for  one  step 

c 

o 

constant 

c 

vector  of  initial  conditions y (0) 

E 

activation  energy  in  the  Arrhenius 
expression 

F 

feed  rate  of  reactant 

F 

-a 

block  entry  in  the  fundamental  matrix 
contributed  by  the  parameter  vector 

F  . 

-xu 

block  entry  in  the  fundamental  matrix 
contributed  by  the  jth  state  vector  s^ 

f 

fractional  conversion 

f  .  (x  ;  a,  t)  • 
3  ~3 

rate  of  change  of  x.  with  respect  to 
the  independent  variable 

Ga  . 

— 3 

block  entries  in  the  fundamental  matrix 
contributed  by  coupling  of  the  state 
vector  Xj  with  the  parameter  vector 

£(y,  a,  t) 

rate  of  change  of  the  total  state  vector 
with  respect  to  the  independent 
variable 

H+ 

concentration  of  hydrogen  ion 

J 

Jacobian  of  cj,  (Z'  a >  t)  with  respect  to 

Z 

. 
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J  . 

-a: 

Jacobian  of  f . (x . ,  a. 
to  a  “=> 

t) 

with 

respect 

J  . 

-XU 

— 

Jacobian  of  f.(x.,  a. 
to  x  “3  -D 

~3 

t) 

with 

respect 

K 

e 

- 

eauilibrium  constant 

L 

— 

the  distance  from  the 
tubular  reactor 

entrance 

:  of  a 

L 

m 

- 

the  total  length  of  a 

tubular 

reactor 

m 

N 

N 

N 

N 

N. 

N 

N 

i 

n 

P 


Cl, 


B 


HCl 


h2o 


Pe 


0 


the  superscript  referring  to  the  cruasi- 
linearization  iteration  count 

the  number  of  moles  of  chlorine  initially 
present  in  the  reactor 

the  total  number  of  moles  initially 
present  in  the  reactor 

the  number  of  moles  of  hydrogen  chloride 
initially  present  in  the  reactor 

the  number  of  moles  of  water  initially 
present  in  the  reactor 

the  number  of  moles  of  inerts  initially 
present  in  the  reactor 

Peclet  number 

feed  rate  of  inerts 

the  number  of  experiments 

the  total  pressure  of  a  point  in  the 
tubular  reactor 


- 

the 

partial 

—  2 

P 

— 

the 

exit  or 

e 

PHC1 

— 

the 

partial 

P 

H2° 

- 

the 

partial 

pressure  of  chlorine 
final  pressure 

pressure  of  hydrogen  chloride 
pressure  of  water 


. 

>  . 
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R' 

A 

R 


r 

S 

max 


s 


v . 

-D 

W 


W' 


the  inlet  of  initial  pressure 

the  partial  pressure  of  oxygen 

the  total  pressure  of  the  reactor 

the  order  state  variable  vector  x. 

( j  =  1 ,  2  ,  .  .  .  n) 

the  order  of  the  parameter  vector  a 

the  ideal  gas  constant  (Kcal ) / (g-mole ) 
(  K) 

2 

the  ideal  gas  constant  (cm  ) (mmHg)/ 
(g-mole) (°K) 

the  relative  step  size 

the  rate  of  reaction  of  oxygen 

the  order  of  the  vector  y_ 

the  initial  maximum  step  size 

the  cross  section  are$  of  a  tubular 
reactor 

temperature 

the  independent  variable  (or  time) 

the  kth  discrete  value  of  t 

the  kth  discrete  value  of  t  for  which 
data  is  given  for  the  jth  experiment 

a  particular  solution  to  the  system  of 
first  order  linear  differential  equa¬ 
tions 

an  intermediate  solution  required  for 
the  jth  experiment  when  the  initial 
conditions  are  assumed  exact 

the  concentration  of  water 

the  weight  of  catalyst  present  in  the 
reactor 


Ill 


w 


X 


ijk 


X 


* 

ijk 


Y 

Z 


Y 


ik 


Z 

z 


a 

Y 

6 


e 


ijk 


a  solution  to  the  linear  systems 
associated  to  the  maximum  operation 

the  predicted  value  of  the  ith  state 
variable  for  the  jth  experiment  for  the 
kth  discrete  value  of  t  for  which  data 
are  taken 

the  experimental  value  of  x  — ^ 

the  pseudo-experimental  value  of  x. 
derived  from  data  perturbation  ^ 

the  state  vector  for  the  jth  experiment 

the  fundamental  matrix 

the  state  vector  which  includes  the 
states  for  all  experiments  and  the 
parameters 

the  predicted  value  of  the  ith  element 
of  y_  for  the  kth  discrete  value  of  t 
for  which  data  is  taken 

the  experimental  value  of  y^ 

the  pseudo-experimental  value  of 
produced  by  data  perturbation 

objective  function 

dummy  vector  required  for  the  maximum 
operation 

the  order  of  the  reaction 

perturbation  constant  for  the  initial 
conditions 

a  vector  which  must  be  greater  than  or 
equal  to  zero  for  all  t  if  the  posi¬ 
tivity  requirement  is  to  be  satisfied 

the  vector  of  errors 

the  deviation  of  x.  from  x. 

13k  13k 
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APPENDIX  A 

Documentation  of  the  Basic  Program  and  the 

Nonlinear  Example 


The  quasi linearization  and  data  perturbation 
procedure  was  programmed  in  Fortran  IV  to  solve  the 
problem  of  finding  the  characteristic  parameter 
vector  a  of  a  system  described  by  a  system  of  differ¬ 
ential  equations  of  the  form: 


dx 

dt 


=  g  (x  ,  a,  t) 


(A-l) 


It  was  assumed  that  n  different  experiments  were  per¬ 
formed  with  the  known  initial  conditions  x_.  (0)  j  =  1, 

2,  ...  n  and  that  some  or  all  of  the  state  variables 

were  measured  at  several  discrete  times.  Given  the 
model  and  the  data,  the  parameter  vector  is  found  such 
that  the  unweighted  least  squares  are  minimized. 

The  Fortran  IV  program,  listed  on  the  following 
pages,  is  broken  up  into  a  main  line  and  five  subroutines 
The  purpose  of  each  subroutine  is  stated  in  the  comment 
cards  at  its  beginning.  The  main  line  program  reads  in 
the  required  data  and  controls  the  data  perturbation 
stepping  procedure.  Comment  cards  in  the  main  line  pro¬ 
gram  defined  each  variable  which  is  read  in  as  data.  All 
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of  the  subroutines  are  general,  except  for  the  sub¬ 
routine  SLOPE  which  defines  the  functions  and  Jacobians 
required  for  the  particular  problem.  The  subroutine 
SLOPE  must  be  supplied.  A  well  documented  example  is 
listed  with  the  rest  of  the  program. 

The  SLOPE  subroutine  listed  is  the  one  required 
to  solve  the  nonlinear  example  which  was  described  by 
the  differential  equations: 


dx^ 

dt 


-ai(xi 


x i  -k  0 

1  e2 


u-x1-x2r| 


dx2 

dt 


,  2 
al(xl 


(l-Xi-x2 


>] 


(A- 2) 

The  input  data  required  by  the  program  may  be 
broken  into  three  blocks.  With  the  first  three  read 
statements,  the  variables  which  control  the  program 
are  defined.  In  the  second  block  the  data  for  each 
different  experiment  is  read  into  memory.  The  third 
block  of  data  is  that  which  is  read  on  the  first  entry 
into  the  subroutine  SLOPE.  This  third  block  consists  of 
the  extra  data  reaired  to  define  the  functions  and 
Jacobians  for  each  experiment.  In  many  cases,  such  as 


, 
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with  this  example ,  the  third  data  block  is  not  needed. 

A  complete  description  of  each  variable  read  from  the 
data  cards  is  given  on  the  comment  cards  contained  in 
the  Fortran  IV  listing. 

The  data  was  obtained  for  this  example  by  using 
equation  (A-2)  with  the  values  of 


II 

1 — 1 
<3 

2 . 0 

k  , 
el 

=  1.8 

a2 

3.5 

k  0 
e2 

=  3.0 

P 

u> 

II 

5.0 

k  _ 
e3 

=  1.0 

for  the  parameters  and  the  equilibrium  constants.  Three 
different  initial  conditions  were  used. 


1" 

0 

0 

(0)  = 

0 

;  x2(0)  = 

1 

;  x3 (0)  = 

0 

(A-4 ) 

The  data  generated  for  each  of  these  initial  condition^- 
is  correct  to  within  ±  0.00005.  It  would,  therefore,  be 
expected  that  the  parameters  extracted  from  the  data 
would  be  almost  identical  with  those  given  in  (A-3) . 

This  was  verified  with  the  results  obtained  from  the 
first  and  third  sets  of  control  data  cards  shown  immedi¬ 
ately  following  the  program  listing  and  input  data. 
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L 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 


c 

£ 


c 

u 

c 

c 

c 

c 


THt  PURPOSE  UF  THIS  PROGRAM  IS  TO  IDENTIFY  CERTAIN 
UNKNG a i\  PARAMETERS  IN  A  SYSTEM  GF  UKL/INakY 
DIFFERENTIAE  EQUATIONS,  WHICH  MODEL  A  GIVEN 
PHYSICAL  b  Y  b  T  E  M  FRUM  WHICH  EXPERIMENTAL  LATA  HAS 
BEEN  OBTAINED, 


THt  PRUBLcM  Is  TKEATEL  AS  A  NONLINEAR  MULTIPOINT 
BOUNDARY  VALUE  PROBLEM  ANu  SUlVLL  USING  A 
COMBINATION  OF  QU AS IL I NEAR I L AT  ION  aND  BOUNDARY 
VALUE  PERT  UbAT I  ON . 


kEAl  U6 ( 10 , 10) , C( 11, 10 ) , XG<30, 10) , TCI  30 , 15) , X8(  30,  15,10 
1  )  ,H,XBAR (30 

1«  15  #10)  »  Xbl ( 30,  15 »  10)f  RPtRAT  IG»RMAXI  , COI 10 ) , X< 10 , 10)  »CH 
1MAX, 

1P(1Q) tT. SUMA  » SOMB 


I N T  E  G  E  R  iVE,NR,RMAX,NUSLL» I C  0 T  M  » 1 0  U  T 1 , NG,NPTCT»NP(30)  , N  P 
II J (30,1 5),  N 

1C  (  b'v  ,15)  » l\  P 1  T  0  I  (30  )  ,  R  S  1  (  3  C  ,15,10)  ,  M 1  ,  M  2  ,  N  C  J  » J  ,  I  ,  II  » J  J  » t\ 
1 , S , LL 

I »  LR  J  M , CK  J i ,  I ChECK ,LKJS,RMAXS, ICUTS 

_ LUMiMlN  Z  A d b/ XC  «  1  ,  n Db  /  Ao3  /  C  ,  N  b  ,  NK  ,  NO  /  A  B4/  UB , X , P  ,  H  ,  T  ,NCJ  , 

ILL/ AblO/TC, 

lXb, ERROR, CHM AX, R MAXI , KS i , NPI J , NC , NP , NP I T 01 , RMAX , NP TOT , R 
I,J,Mi.»M2»Ii 

1 ,  J J  i  S ,  ICHEuK/AB4/LNJM, lKJ I /A 66/ CU 


TO  USE  THIS  PROGRAM  THE  DIFFERENTIAL  EQUATIONS  MUST 
oE  NORMALIZED  bUCH  THAT; 

1  nCNl  OF  THE  STATE  VAK I A BE Eb  EXCtLDS  10. 

2  THE  INDEPENDENT  VARIABLE  V  ARK  I Eb  FROM  ZERO 


C 

t 

c 

C 

C 

C 


TU  A  VALUt  NuT  EXCEEDING  100,  WITH  THt  INITIAL 
C0NU1T10N  BE  I  No  blVEN  AT  ZERO. 

3  THt  ELEMENTS'  UF  THE  PAR AMt Tt  R  VECTOR  ARE  UF 
THt  bAMt  APPROXIMATE  ORDER  OF  MAGNITUDE. 


C 

L 

C 

C 

c 

> 

c 

c 

c 


THE  MAIN  LINE  PKOCAM  mAS  THt  FUNC5 IONS  OF  kEaDING 
IN  1  HE  DATA  SETS  AND  CONTROLLING  THE  PtRTURBtD  DATA 
STEPPING  PROCEEDURE.  THE  DATA  A*t  WRITTEN  EXACTLY 
mj  THEY  AKt  READ.  *  BEGIN  PHASE  I  *  Is  THEN  WR I TTEN 
CUT  SIGNIFYING  THE  BEGINING  OF  TmE  STEPPING 
TROCEEUUKE.  II  IS  ALSO  WRITTEN  GuT  EACH  TIME  THE 
bTEF  SiZt  lb  LclKlasED  AND  THE  PkUCtEDUKE  lb 
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RE~  I  fSi  i  T I  A 1  ED.  A3  THE  BEGIN! NG  OF  EACH  S3EP  THE  ST  tP 
i\c.  IS  WRITTEN  ALUNU  WITH  THE  CURRENT  ESTIMATES  OF 
THE  PARAMETERS  AND  I h E  SUM  LF  THt  E  k  u  r<.  b  S^uAkEl). 
'oEGIN  PHASE  II*  IS  WR 1 1  TEN  WFIEN  THE  ACTUAL  UATA 

-  POINT  S_  ARE_JiS.EC.  PHAS  E  II _ IS  TERMINATED  WHEW 

llNVERGENCE  IS  OBTAINED  UR  WHEN  TEN  SUCCESSIVE 
i  3  ErAT  IuNS  FiAVE  btth  COMP LE TED  ^Ilriuur  CONVERGENCE 
CR  WHEN  PHASE  I  HAS  BEtiM  RE-ENTEkEl)  BECAUSE  OF  AN 
ExCESIVELY  LARGE  CHANGE  IN  THE  PARAMETERS  IN 
SuCCESbiVE  ITERATIONS.  IN  THE  LATER  CASE  'BEGIN 
_ _ PHASE _ I*  IS  AGAIN  WRITTEN  UUT. _ 


auxiliary  scbkuui ines  rewuireu  a  r  e  j 

1  GENER  —  UuEl  TO  GENERATE  THE  SOLUTIONS  FOR 
APPROXIMATE  PARAMETERS  NEEDED  FOR  THE  STE PINb 
Pr\L  l  1 1  uUR  c  AND  FOR  GENERATING  ThE  FINAL  SOLUTION 
SC  THAI  THE  SUM  LF  THE  ERROkS  SlUaREU  MAY  BE 
CALCULATES  . 

2  uUASI  —  USED  10  CONTROL  THE  GU AS  I L I NE AR I L AT  I ON 
ITERATIVE  PRQCEEOURE. 

3  1NTE  -  CALLED  FROM  QUASI  AND  USED  TO  INTEGRATE 
THE  SYSTEMS  OF  DIFFERENT  I AL  EQUAT  IONS  ENCOUNTERED 
a i IH  THE  QUASI LINEAR I Z ATI  ON  PRQCEEDURE. 

4  LEAST _ ~~  CALLED  FR  QM  ^UAbl  AND  UbED  Tu  PET  ERMINE 

THE  LEAbT  bbUARE  MATRIX  EQUATION  WHICH  IS  SOLVED 
IN  THE  GAGSb  ELIMINATION  SUdkuUTINE. 

5  b A b S b  -  CAlLcu  FROM  LEAST  AND  Io  UbED  AS 
DESbRiBLD  AbOVE. 


b  slope  -  Called  from  uEner  and  i n t e  and  is  used  to 

EVALUATE  THL  FuNCTIUNS  AND  JAlUbIANS.  THIb 
SUBrUUTINE  FiUbl  Be  bUPPLiED  AlUNU  WITH  THE  DATA. 

THc  FURMAT  OF  THIS  SUbRUUT  iUE  IS  EXPLAINED  IN  THE 
SAMPLE  uIVEN  IN  THIb  LISTING. 


THc  FOLLOWING  I NT EGeRS  AKt  READ  IN  ON  A  I6lb  FORMAT. 

NE  THE  NO.  UF  S3  ATE  VARIABLcS.  NE<11 

NK  THE  NO.  OF  PAKAMATeRS.  NK<11 
C 

L  R M A X  THE  MAXIMUM  NO.  CF  COMPLETE  ITERATIONS  TO  BE 

C  MADE  A 1  EACH  STEP.  RMAX<jLl 
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NDS  THt  NG.  GF  DATA  SETS.  NDSOl 

LL  PRINT  GUT  CONTROL  FOR  STATE  VARiABLE  PkUFILES; 
_£UK_  tXATiPLJ^J  E  LL  =  2,  THEN  EVERY  SECOND 
INTEGRATION  STEP  IS  PRINTED.  A  *  IS  PRINTED  TO 
THE  LEFT  GF  EVERY  VALUE  GF  THE  INDEPENDENT 

vakiadle  fgk  which  data  are  given. 

I  Cu  TM  THE  PAX  I  RUM  NU .  GF  STEPS. 


LKJM  PRINT  SGPRESS1CN  CGNTRGL. 

LK J  M  =  0  PRIMING  GF  STATE  VARIABLE  PROFILES  IS 
DEPRESSED. 

ER J  M  =  K  STATE  VARIABLE  PROFILES  ARE  PRINTED 
IN  PHASE  I  AND  II  FUR  EVERY  K  TH  ITERATION. 

L  K  J  R  - — K  STAT  E _ VARIABLE  PROFILE  S  a  R  E  PR  IN  T  c.  U 

FOR  EVERY  A  TH  ITERATION  IN  PHASE  II  ONLY. 


RE Al 1 5  » i  )  N E  ,  NK , RM  A  X , ND S  9  LL ,  ICOTM»LKJM 
wk I T  E ( 6 , 1  i  NE  ,NK , RM AX , NUS , LL , I CUTM, LKJ  M 
LKJ.S  =  LMUM 

RMAXS=RMAX 
I  CU  T  S  =  IGLTM 
Nl=I 


CU _ THE  IN I T IAL  GULSS  FuK  THt  PARAMETER  vECTuR  IS 

READ  IN  UN  A  B  E 1 6  •  o  FORMAT. 


REAL (5,11)  ( CGI  I ii , I  1=1, NK) 
WRITER,  11) ( Cu(  I  I)  ,  I  1=1,  NK) 


THE  FOLLOWING  CONTROL  VARIABLES  ARE  READ  IN  UN  A 
4  E 1 6  •  6  FORMAT. 

ERROR  f HE  MAXIMUM  aMQUNI  OF  RELATIVE  EKROR  AEEGwEO 
_ I_N  _  TH  E  DESI  RED  P  A  R  AML  T  E  R  S  . 


RMAXI  THE  INITIAL  SlEP  SIDE. 

OHM AX  THE  MAXIMUM  ALLOWABLE  RtLATlVE  UHANuE  IN  THE 
PARAMETERS  FuR  UNc  ITERATION.  IF  ChMAX  IS 

_ EXCEEDED  kMAXI  IS  HALVED  AND  Pn«oE  I  IS 

RE-INI 1 1 AT  ED* 


lhmin  the  minimum  allowable  relative  change  in  the 
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PARAMETERS  FUR  ONE  STEP  BEPORc  RMAX  I  IS 
TRIPLED. 


READ (5  ,  11J  ERRuk i KMAXI , CHMAX , OHM  IN 
WRITE! 6, 1 1J  ERROR, RMAX  I ,CHMAX ,CHMIN 
II  FORMAT !  5E  16. 6 ) 
l\PTGT  =  C 


_ THE _ DATA  FUR _  EALh  EXPERIMENT _ ARE  -READ.  I N _ 

SEQUENT  I  ALLY  AS  GIVEN  BELOW.  ALL  INTLuEKS  ARE  READ 
IN  ON  A  1615  FORMAT  AND  ALL  OTHER  DATA  ON  A  8FI0.6 
FORMAT.  THl  DATA  1 5  PRINTED  OUT  EXACTLY  AS  IT  IS 
R  c AC  IN. 

_ I _ LS _ THE _ LnELE-P-E ND ENT _ VARI ABLE. _ 


DO  101  1=1 »NCS 


XO(I»J)  TtiE  INT1AL  CONDITION  FOR  THE  J  TH  STATE 
_ VAR  I  ABLE  FOR _ THE  I  Th _ EXP  Ek i ME  NT  . _ 

READ!d,2I  (XCllt J) » J=1,NEJ 
WRITHE*  2J  (  XL!  I  ,  J)  » J=  1  *  NE  i 

NP(1)  TH t  NO.  OF  VALUES  OF  I  FOR  WHICH  DATA  AkE 
_ GIVEN  IN  THE  I  TH  £ XPERT  MENT . 


REAL (5,11  N  P ! I  ) 
aRI TEib* 1 JNP! I ) 
PI  =  N  P  t  I  J 


KFIjC  I  * J)  THE  NO.  OF  DATA  Pul  NTS  WHICH  ARt  GIVEN 

FOr<.  ThE  J  TH  VALUE  OF  T  FOk  wHiCH  DATA  ARE 
GIVEN  FOR  IhE  1  TH  EXPERIMENT. 


R  E  AO  l  5  *  1  i  (NPIJU,JJ,J=1,M1) 
WR  I  T  t ( 6  » 1 ) (NPlJ(ItJ)  *  J  =  1 1 Ml ) 


Tuli»OJ  TFE  VALUE  OF  THE  J  IH  VALUE  OF  T  FUR  I G H 
DATA  ARE  GIVEN  FOR  THE  1  Th  EXPERIMENT. 

Rc  A  D ! 5  >  2  )  (TCtl*J)  »J=1*M1) 

WR ITE(6*2)!TC!I*u)  »  J= 1  *  M 1 J 

NCUiJJ  THE  Nu.  OF  INTEGRATION  POINIg  REQUIRED 
BETWEEN  1 C (  I  *  J  )  AND  TG(I~1»J)  aITH  T  C ! 0  »  J ) 
BEING  TAKEN  AS  aERO. 
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C 

ROAD ( b , 1 )  ,  J  =  1 ,  M.  1 J 

to  R I  T  E  {  6  »  1)  (  N  C  (  i  f  J  )  »  J  - 1  » M  1  ) 

INP1  TCT  t  I  )  =0 

_ QQ  1QQ  J-1.M1 _ _  _ 

M2«NPI J ( I , j ) 

J\  P  l  T  G  T  (  I  )=NP  l  TOTi  I  )+M2 
C 

C  Ks  I  l  I  y  J  r  K  )  I  HE  sTaTE  VARIABLE  INDEX  OF  1  HE  K.  TH 

C  DATA  POINT  AT  TO <1,  J). 

C _ _  _ 

REACH  5, li  IKS  I  (I , J  ,Ki ,K=1, M2) 
WKITE(6flJ(KSnifJfK),K=ltM2) 

C 

C  Xdi(I»J»K)  1  HE  VALUE  LP  THE  K.  TH  DATA  POINT  AT 

0  TO (  I »  J )  wITH  THE  INDEX  KSI(I»J,K). 

jC _ _ _ _  _ 

READ  15*2)  (Xfil(i,J,K),K=l,M2) 

100  wR I  Ic i c , 2J( Xdl  (  I , J  f K J , K=1 ,  M2 } 

10  1  NPTOT-=NPTOT  +  NPI  T0T  1  I  ) 

C 

c 

c _ XfaE  TOTAL  NO.  OF  DATA  POINTS  IS  WRI  TEN  OUT . 

C 

C 

wR I T  E ( 6  *  22 )  N PICT 

22  FQRMATl'QTHE  TuTAL  NO.  OF  DATA  PLINTs  IS  »,I5) 

1  FORMAT i 161 IX , 14) ) 

2  FORMAT <611X»fG.6) i 

C 

C 

C  IF  THE  Total  NO.  OF  Data  POINTS  EXCEEDS  5GQ  THEN 

0  THE  PROGRAM  IS  TERMINATED. 

C 


IF (NPTCT.LE.500J  GOTO  120 
WRIT£(6»21) 

21  FORMAT  (  1 ONO .  OF  DATA  POINTS  EXCEEDS  DOOM 
STOP 


0  PHASE  I  IS  INITIATED. 

c 

c 

120  TOUT  I  =  C 
114  1CHECK=1 

_ WR  1  TE  (  6  f  18  J 

18  FORMAT ( *  OB EG  IN  PHASE  l*) 
LK jM=LK JS 
RMAX-KMAXS 
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ICGTM=ICCTS 

C 

c 

C  GENERATE  THE  SOLUTION  FROM  THE  APPROXIMATE 

_ C _ PARAMETERS. 

C 

C 

110  R  P  =  0 . 0 

DU  107  I  —  1 ,  !\  K 
10  7  C( 1 , I ) =Co (  I  ) 

_ SUMB-0 . 0 _ 

DC  103  1=1, NOS 
T=0 .0 

DU  104  1  I  =  1 ,  N E 
i  0  4  XI 1  ,  i  i  )  =  X G ( 1,11) 

M 1  =  N  P (I) 

_ DO  1C3  J=1,M 1 

l\C J=NC { I , J ) 

H— (  rci 1  ,  Jj-T )/NCJ 
GALL  GENER 
M2=NPi J l I , J) 

CO  103  K=1,M2 
_ S—  R.3  1 1 1  f  J  ,  N ) 

XBAR  {  i ,J,KJ=X(1,5) 

SUM A= A ESI XB1 I i , J ,K )-XBAR l I , J, K) ) 

SU MB = S CM b+SUMA*#2 
1 C  3  IF(SCMA.GT.RP)  RP=SUMA 
ICGTI-ICUI 1  +  1 

_ c _ 

c 

C  STAkT  STEP  ICUfi  AND  CHECK  TU  SEE  IF  XCGTM  HAS  BEEN 

C  EXCEEDED. 

c 

c 

WRITEC6. 14)  ICUTI 

14  FORMAT  I  * C* , 10X, *  BEGIN  STEP’,15) 

KRITE(6»23i  (I,CC(1),I— 1,NK) 
vyR11E(c,19)  s  g  mb 
1F( ICUTI .LT .  ICUTM)  OgTU  119 
nR I T  E ( t , 20 ) 

2 0  FORMAT  t  *QST£P  COUNT  EXCEEDED*  ) 

STOP 

119  RAI 1G=1 .C-RMAXi/KP 
C 
C 

C  CHECK  10  DETERMINE  IF  PHASE  II  MAY  be  ENTEkcU. 

C _ 

c 

IF(RP.GT.RMAXI)  uOTU  111 
in  HI  TEIb, 17) 
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17  F  t_RM  AT  (  *  lutOlN  PHASE  I  1  ■  i 
R  Ml  A  X  =  1 0 
R A  T I G  =  G  .  C 
LKJM-IAdSlLKJM) 

i£6..rM.=-  ilgtm 

c 

c 

C  UcTERMiiNc  P  S  E  U  U  C  l>  A  T  A  S  £  T  5  • 

c 

c 

_ 111  DO  106  1  =  1, NLS 

P1  =  NPU  ) 

CO  106  J= 1 , M  1 
M  2  =  IN  P  I  J  (  I  ,  J  ) 

OG  106  K=  1 ,  M  2 

XU  IfJ»KJ  =  XblC  1  ,  J  ,  K)  +RAT  Iu*i  XbARi  I  ,J,K)-XBlUi'J,K)  j 

_ 106  CONTINUE 

C 

C 

C  EiN  F  ER  *UASlLlINiERIZAT  I  GN  ITERATIVE  PRuCEEOUKE. 

C 

C 

D  A  Li  w  G  A  S  i 

C 

G  CHtGK  FUR  PHASE  I  RE-ENTRY. 

C 

c 

IF (  IgHECK.lt .0)  GCTC  1 14 

c _ 

c 

C  GHEGK  PGR  PRCGKAR  COMPLETION. 

C 

C 

IFUCUTN.lt  .0)  GET  G  112 
SUM A  =  C  .  0 

CO  lib  1  =  1,  INK 

118  SUMA  =  SUMA  +  ABSUC(  1  ,  I  )-C0  4I  J  )/C  U  ,  I  J) 

C 

G 

c  check  tu  determine  if  step  sizl  should  be  increased. 

c _ 

c 

IF ( SuMA.LT .CHM  IN)  RMAX1=3.0*RMAXI 
GUT  C  ilO 
C 
C 

C_  THE  FINAL  PARAMETERS  ARE  USED  TO  DETERMINE  THE  SUM 

C  LF  THE  ERRORS  SUUARED  AND  ThE  SGLUTIUN  IS  WRITTEN. 

C 
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*  1  =  JL  *  J  J 
^  ■  IUt*i.  C  i  i 

•J 

i  •  •  j  Ju  >-  J  -  i  :.  J  I  „  1  l  J  I  J  .,1  _  I  .U  jfiJ 

-  I  *  i  i  l  ,  I  J  .  I  .  c.  J  11 
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112  DO  115  1=1, NK 
115  C( 1. 1 ) =CC  (  i  ) 

SUMB-0 .0 

CL!  llo  1  =  1, MCS 

T=Q  »0 

CO  117  11  =  1,  fvE 
117  X<  1  ,  II  )  =  XU  I ,  II  J 
M  =  NH  m 
CC  116  J  =  1 »  M 1 
N  C  J  =  N  C  (  I  ,J) 

H=  IT L (  1 ,  JJ  -  I  J/r.LJ 
CALL  0 EN E K 
M2=NPI  J(I  ,  j) 

CO  116  K=  1  ,  M  2 
S=KSI(l,b»K) 

llo  SUMB=SUMB+  (XB1(  1 ,  J  ,K)-X<  1  ,S)  )**2 
IF(  1 CH  ECK  « E  U  .0  )  kRITE( 6,23) 

23  FORMAT  (  *  OCCNV  ERGEN  CE  TO  THE  uESIREC  ACCURACY  AAS  HOT  *  , 
1GBTAINEO  AFTER  1C  ITERATIONS  IN  PHASE  II1) 

v\  K  1  T  E  i  o  »  2  4  ) 

24  FORMAT  I *CTnE  FINAL  SOLUTION  OBTAINED* ) 

WKIfEi6,25J  11,0011)  ,1  =  1  , N K ) 

25  FORMAT  (  *  Q  _  K  *  ,  I  L  »  *  ..=  •  ,E16.6) 

write(6,19)  some 

19  FORMAT C C* ,9X  ,* THE  SOM  OF  THE  SQUARES  OF  THE  ERRUkS  IS' 

1 »  E 14 . 4  ) 

io  FORMAT ( 1 C  AFTER*, 15,*  STEPS') 

STOP 

_ END _ 
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SubRGuJINE  w  0  A  i>  I 

il^tbLK  KSI  (30,15,  10)  »NPIJ(30.  15)  ,NC(30,15)  ,NP(  30)  ,NPI  T 

i or (30) ,ne, 

1NK»RMAX»NDS»NPTCI  ,  R  »  I»J,PlfN2,  INDEX,  1CGUNT»II,JJ,S»LL,N 

LC 

1  ,  L  K  J  iv  ,  Ll\  J  I  ,  IDHECK 

REAL  UBt  10, 10) ,Xl 10, 10)  ,P( 10) ,TC( 30,151 ,XO(30, 10) ,XB(30 
i  ,  1  b  ,  1 Q  )  ,  T  ,  H 
1 ,Ci  1 1 , 10),CR ( 10 ) 

REAL*d  B { 500 ,  10 ) , Q ( 500  ) 

_ C  LiM  R  ON  /  Abb/  CR  /  ABI/b  « U  /  Add/  XG  »  1  «  NO  j/  Ab3/C  tNfctNi^iNU/  Ad  4 

1/UB , X , P j  H i 

IT  ,  N  C  J  ,  L  L  /  AblO/ICyXB,  ERRGR,CHMAX»RMAXI  »KSl,NPIJ»NC»NP»NP 
liTLT ,RRAX , 

1HPT0  T  y  Rf  J',  Ml  ,  H2:,  II,  JJ,S,  I  CmeCK/ AB4/i.KJH,LKJI 
C 
L 

C  THE  PURPOSE  OR  THIS  3UBRG0T INE  IS  TU  CONTROL  THt 

C  QUASILINEARIZAT  ION  ITcRaT  IVt;  PROCEEDURE. 

C 

C 

L  K  J  I  =  0 

CO  100  K  - 1 ,  mvAX 

L(\di-LKJ1  +  1 
I  N  D  E  X  =  C 

DO  104  IsUMS 
S  =  R-  I 

IHLKJi  •NE.LKJM)  GOTO  113 

frjRI TE ( b , 5  J  S  ,  I  ,  NP  I  TOT  (  I  ) _ 

DO  1 12  11  =  1, Nb 

112  WR I T  E ( 6 , o )  I  I r  X C (  I ,  I  I  ) 

to R I T E ( 6 , 8 1 )  ( C { K , I  I )  , 1 1  =  1  , nK } 

WR1  TE( 6,  7- ) 

b  RURiMAl  (  1  5  H 1 E  0  R  I TERA 1  I  ON  , 14,11  H  DATA  SET  ,14,*  WITH  1 

_ 1,14,*  HJIMS  1  /  1  U  f  hE  INITIAL  ClNDiTIuN  IS1/) 

0  FORRAT(10X»lHX,I2f2H  =  , F 1 0 • 6 ) 

71FGRMAT ( 1HC, 14X , 1HT  f8X,2HXl ,8X, 2HX2,8X,2HX3 ,8X, 2HX4,8X, 
12hX5,8X,  2HX6, 8X,  2HX7,8X,  2HX8  ,8X,.2HX9,8X  ,3H.X10  /  ) 

8  FORMAT 419HOTHE  CONSTANTS  ARE  /10X,8E16.6) 

C 

C _ 

0  THE  INITIAL  CONDITIONS  FOR  THE  i  TH  DATA  SET  ARE 

C  INITIALIZED  AND  T  IS  SET  TO  ZERO. 

C 

C 

113  T=0  •  0 

_ DO  105  11=1, NE _ 

DC  1  Go  J  J—  1  , NK 
lCb  UB( II ,JJ)=0.0 
DO  10  7  J  J=  1 »  R 
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10  7 
10b 

X( JO  ,  I  I  )  =XO (  1,11) 

Pill  )  =  XC4I, Ili 

M 1  =  N  P (I) 

DO  13A  J  =  1 ,  M  1 

NC  J  =  Nl  I  1  ,  J ) 

C 

c 

H=ITC( I, J ) -T ) /FLOAT (NCJ) 

c 

L 

C 

ThE  SObROoTINE  INTe  IS  CALLED 
MATRIX  DIFFERENTIAL  EQUATIONS 
UPON  RE  TURN  T=  TC { I , J )  . 

10  INTEGRATE  The 

FROM  T  TU  TCI  I , J )  , 

c 

c 

CALl  INTEiR) 

9 

IFILKJI.NE.LKjM  llTO  IIh 
WRITE(6,9)  T  ,  1  X ( k  ,  1 1  ) ,  I  i  =  l  ,NE ) 
FORMAT ( 7 X  ♦  1H*  .  I0F 10. b) 

c 

c 

L 

SET  UP  ThE  LINEAR  EG AT  IONS  AND 

CALL  SUBROUTINE  LEAST 

TU  DETERMINE  THE  IMPROVED  ESTIMATE  OF  T Hl  PARAMETER 
VECTOR. 


1 1 A  M2=NP I  0 (  i  ,  J) 

Cl  110  K=i,K2 
IClU\T-=K  +  INDEX 
3—  KSi  (  I  »  J  »  K  ) 

_ C  I 1  l  C  U  N  T )  =  X  ti  1 i  ,  J,k)-P( S) 

DO  110  J  J  =  1 »  NK 
i  10  E(  I  COUNT  »vJj)  -Ub  I  S  ,  J  J  ) 

10h  INCEX=  INDtX  +  P2 

IF( LKJM. EC.LKJI )  LKJI=0 
CAlL  lEAST(NK,NP10T) 

Sum a=q  .  o _ _ 

DU  ICS  11=1, NK 

SUM  A=s  UM  A+ A b  S I C ( R ,  1 1 ) -CR (II) ) / ( ABS l C ( R ,  I  i)  ) +0 . 1 ) 
ICS  Ll R  +  i  ,  1  I  )  =CK  1  1  I  ) 

IF1SUMA.CJ.lE MAX)  CLTC  11b 
I F I  SUM A • L T . ERROR )  GOTO  lib 

IF  (j%  .El*  1 0  )  ICHECK^C  _ _ _ 

1FIR.EQ.10)  GOTO  lib 
103  CONTINUE 


THE  PARAMETERS  ARE  WRITTEN  OUT  AT  ThE  END  OF  THE 
S  1  l  P  IF  k  M  A  X  IS  EXCEEDED.  T  ri  I  j  OCCURS  w  H  E  N  THE  ERROR 
SPEC  I F 1C AT IUN  IS  NOT  MET  IN  RMAX  ITERATIONS.  THIS  IS 
C  NUT  TOO  CRITICAL  IN  MOST  CASES  AS  ACCUKATt 

C  SOLUTIONS  TO  EACH  STEP  AkE  l\GT  NECESSARY. 
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.  t  , 


l  *  *  v  V  .  -  I  .J  i  )  \  l  i  ,  i  ■ 


lot  1  i  <  -i  “  1 

t  .  -  /t 

■  ..  i,  1  i  . 

i  t  L>  t  •  )  *  &  1  * 

i  ,  I  ~  ,  t  tit  -  t  .  L-  i  I  J 

•  -i  •  ^  i  J  j.  i 

-  ,  L  J  t  <  L,/>  i  )  1 1 

.  H  t  1  -  l  i  i-  .. 

t  t  t  .  J  / 

i  .  /  -  l  l  i  *  I 

\  .  I  .  'il  ii 

.  .  u  i  i  I 

i  ...  .  J  .  i 

.....  I  i  i 

J  /i  l  I  J  t  . 


ill  v  t  -»  .  ; 

i  .  J  c  ..  I  .  , 
l  ■  .  I  A  1  1  l 

« 


R=RF AX+ 1 

RRITn(fc»lC)  RMAX 

toRi  I  E  ( 6  #  8  )  (C{K.II)»Il=ltNK) _ 

10  FORMAT (7HCAFTER  ,lb,Uh  ITERATIONS) 
115  RETURN 
lib  rMAX1  =  RRAXI*C  .  b 
DO  117  11=1, INK 

117  CR  l  1  I  1=0 (  1 ,  I  I  ) 

1CHECK=- 1 

118  RETURN 
END 
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SUBROUTINE  GENER 

REAL  UB<  10,10) ,P( 10)  ,X( 10,10) ,XA< 10,10)  ,KP (4,10)  ,G( 10)  , 
1  T  ,  T  A ,  H 

1,JX< 10,10) ,JK(10, 10) ,C{ 11, 10) 

I  NT  EGER  NCJ,NE, NK ,  LL ,NO,N, II ,  I 

-COMMON  /J&1LL+UE,  NK  t  NO  / A  R4 /U1L,  X  ,  P  ,  H  «J...  Nr.  J .  13  /AB  5  /  .1  x  ,  .IK  , 
1  XA  ,  G 
C 

c 

c  THE  purpose  of  THIS  subroutine  is  tg  generate  THE 

C  SOLUTION  OF  THE  SYSTEM  OF  DIFFERENTIAL  EQUATIONS 

L - GJ-VFN  A  SFT  OF  PARAMETERS.  T hf  4  th  (ORDER  rmndf _ 

C  KUTTA  PROCELDURE  IS  USED. 

C 

C 

DO  100  N  =  1 ,  NC  J 
DO  101  11=1,4 

- GOTO!  in?.lni,i  r>4,  i  ns  )  T  II _ _ _ 

102  TA=T 

DO  106  I = 1 , NE 

106  XA{ 1 , I )=X( 1 , I ) 

GOTO  107 

103  TA=T+Q.5*H 

- DO-  1,0ft  1  =  1 ,  NE _ 

108  XA{ 1 , I )=X( 1 , I ) t0.5*KP(l, I ) 

GOTO  107 

1X14  DO  109  I  =  1  ,  NE 

109  XA ( 1 , I )  =  X < 1 , I ) +0. 5 *KP( 2 , I ) 

GOTO  107 

105  TA=T+H 

DO  110  I  =  1 , N  E 

110  X A {  1 , I ) =  X ( 1 ,  I  ) +KP (  3,1) 

107  CALL  SLOPE (1 , TA,-1 ) 

DO  101  1=1, NE 

101  KP ( II , I )=H*G ( I ) 

_ I^LLA _ _  _ _ 

DO  100  I=1,NE 

100  X(  1  ,  I  )  =  X  (  1  ,  1  >M  KP(  1  ,  I  )  *2.n*(  KP  <  2  ,  I  )+KP<  3,  I  )  )  +KP  {  4,  I  )  )/6 
1.0 

RETU RN 
END 
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bUbkJUT  1  N  E  1  N  T  c  i  R  ) 

R  E  h  L  OB  (  10,  10  )  ,P  (  iU)  ,XQO  ,10  1,  dBa(  10 ,10)  ,PAt  10) , XAI 10, 1 

IC  )  ,KA(4j1C, 

1 1 0  )  ,  KU fi {  4  »  i  0 , 1 0 )  , K P  i 4 , 10  )  ,  G  l  1 Q  )  ,  T  ,  T A  ,  H  ,  S UM A  ,  SU H B  ,  J X(  10  » 

1_101j _ _ 

lbM  L  ,  1CJ  ,  C  (  11 ,  1C) 

INTEGER  NO  J  ,  R  ,  1  ,  J  ,  i  I  ,  NE,NK,L,LLtS,N,NO,LKJP,LKJI 
C  uH  h  v-iM  /  A  b3  /  C  ,  N  t ,  N  R  ,  N  u  /  A  D  4  /  U  B  ,  A ,  P  ,  H »  T  ,  N  C  J ,  LL/ABb/JA,  Ji\, 
1  XA  ,  G 

1/Au4/LKJR,lKjI 


fHh  PURPOSh  UP  THIS  subroutine  is  tu  GENERATE  the 
S...Q.LUJJiM  0£  IU t  SXSJLEM  £LE  MATRIX  DIFFERENT  iai 
equations  required  for  the  SOLUTION  up  THE  linear 

BOUNDARY  VALUE  PROBLtPb  ENCOUNTERED  IN  EACH 
X_F  E  R  A  T  1  u  4 _ Op  THE  wC'Mb  IL  iNtAR  1  ZA  1  1  ON  PRQCEEoUKE  ,  T  HE 

initial  guess  is  albu  generated,  the  4  th  ordlr 

a  RuNGE  KUTTA  PRGOElDORE  lb  USED. 

C 

C 

L  =  0 

_ CO  1 0  Q  _ _ _ _ 

CO  101  11=1,4 
GoTC(  102,103 ,104,105) ,  II 
i  C  2  TA  =  T 

CC  109  1 s 1 , N  E 
PA ( I j  =  P  (  I  ) 

- CU _ 108  J  = 1 , R _ 

10 d  XA ( J , I j~X { J  ,  I ) 

Du  109  J  —  1 ,  i M  K 
1.C  ;  U  &  A  i  i  ,  J  )  *  G  B  l  1  ,  J  ) 

GLTC  110 
1 C  3  TA=T+G ,5*H 

CO  113  1=1, NE _ 

PAi I  )  =  P(  I  )+Q.B*KP( 1,1) 

DO  1  i.2  J  =  1 , R 

1X1  x A ( J  ,  1 J  =  A { J , I i+Q.5*KX( I«J. I ) 

CO  ii  j  J  =  1  ,  K 

113  UB A  { 1 ,j)=Ub(  IfJ)+0.5*KUB(l,l,Ji 

_ GOTO  110 _ _ 

ICh  DC  114  1=1, NE 

PA( I  )  =  P  1  I  )  +0 .5 ^KP (2,1) 

CO  LL6  J  =  1 , R 

116  X  A  l  J  ,  1  )  =  X  (  j  ,  1 )  -t-  G  .  D  X  (  2  ,  J  ,  1  ) 

DO  114  J  = 1 ,  N  K 

11  A  UBA  (  1  ,  j)  =Up  (  I , j J  +Q  .  b*RUB  (  2  ,  I  ,  J  ) _ _  _ 

GOTO  110 
1C  3  TA  =  T  Hi 
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PA l 1  )  =  P (  i)  +KP(3, I ) 

_ _ _  Lb  JUL9  U  =  1  j R 

119  XAiJ,I)=X{J,I }+KX< 3» J 1 1 ) 

OG  i  i  7  J= 1 »  N K 

_ 117  UbAI  I  »  J  )  =  bb  t  i  « J  )  +  Kbb  (  3  »  I  i  J  ) _ 

11C  lalL  SLOPE  { R  ,  T  A  ,  1 ) 

1 F ( R  .  M  E  .  1 )  GbTb  lib 
OG  LL5  l  -  1 1  f .  c 
SUMA=o  (  1  ) 

1 1  ^  KX I 1 i , 1 , I i=h*SUMA 
1 1  o  DG  12C  1=1, Nc 

S  U  M  A  =  G  l  1  i 
LG  12  i  J  =  1 ,  N  K 
SUMB-JKtl ,  J) 

S  y M  A=  S  0  M A- J  K ( 1 , J ) *C ( R , J I 
CG  1  22  $  =  1  »  N  E 

_ 122  SUMb=GUMb+ JX i I  »S)*GbA(S,U) _ 

121  KUB  ( I  1 , i ,  J  )  =  H*SUMB 
LG  i 30  J— 1 , N E 

13b  SGMA=  SUMA+ JX  i  i,J 1*  I  PAL  J  )  -XAlR  ,  J  )  l 
12b  KP ( 1 1 , 1 )=H*SGMA 

IF ( R . EC . 1 )  uCTG  101 
DU.  124..  .S=.2a R 

CALL  SLOPE! S-1,TA»1) 

IF ( S  .Nt • 2 i  GCTU  107 

UG  ill  I  =  i  ,  G  L  — - - 

SUMA=o  I  1  ) 

KXUiiLtl  )  =H*SUMA 

_ ill  bbis.iGbc _ 

iC7  DG  124  1  =  1,  i\  £ 

SUM  a  =G  (  I  ) 

DQ  12.5  ...  J=.1..,..M.E _ _ _ — 

12  3  S  U  M  A  =  S  U  M  A  +  J  X  l  1  »J)*(XA(5,J)~XAlS~l»J)  ) 

DG  12o  J  =  1 » N  K 

12c  SUMA-GuMA-t-JN  (  1  .  J  )  »  i  L  (  S  ,  J  j  -C  1  S~.Li.JJJ _ 

124  KX  (  I  I  ,  G,  1  )  =14*  SOM  A 
101  CGi\  TINGE 

„ _  DO  1  <_  i  i  =1  J  .N..E _ ___ _ _ _ — — — — — - - 

P(l)  =  Pm  +  tKP<l,IJ+2.0JMKPl2,IJ  +i\P  <3,1  )  )  +KP(  4,  I  J  )  /6.0 
Lb  127  J=  1 , NK 

_ Ob  I  I  ,  J  )=Uu  (  1  ,  J  J-KKU3  (  1  .  I  ,  J  )  +2 .0*  (KOb  (2  ,  1 ,  J  )  +  KUb  (  3,  1 ,  J  )  ) 

1 +  KL  b l 4  ,  I »  J  I 
1)/6.G 

b 

b  CtitCK  TG  lETLRMINL  IF  THE  FUNDAMENTAL  MAlRlX  IS 

0 _ b GONG  EL  .  IF  NUT  EXECUTION LG TERM!  .NAI  EU  « LE LHxlKL - 

C  IS  A  FAILURE  AT  THIS  POINT,  IT  lG  lIkElV  THaT  TnE 

C  I  NT  EGR  AT  I  C  IN  PRubELDURt  IG  U  !\  S 1  A  b  L  E  FOR  THIS  SYSTEM 

b  WITH  THE  GIVEN  NQ .  Of  INTEGRATION  POINTS  NC( I, J ) 
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L  Ai^G  MTh  [  H  E  CURRENT  ESTImENT  OF  THE  PARAMETER 

£_ _ _ _ ...VECTOR. 


- If(A8SlUti<I.  JJl.LT.l.OE  50)  GOT  Q  l?7 _ 

RRI  T E ( 6 1 2 ) 

~  FORMA  1  i'CiHE  INTEGRATION  UP  THE  FUNDAMENTAL  MATRIX  iS*» 
i ’  DABO UNPEG* i 
STOP 

12  7  C  G  N  T i N U  E 

_ DO  129  S=  1  . R  _ 

CO  120  I  —  1 f  N £ 

12  U  Xl$ili=X(G,I  i  +  {KX(lfS»li+^.0^(KX(2»S,Ij+KX(3»S»n)+KXC4 

i ,  g  ,  i  u  /  o . : 

1  =  1  A 
L  =  L+  1 

- IFa.Nt.LL)  GOTO  1  00 

L  =  C 

IFIUJi  .NE.LKJMIGCTC  100 
IF(i\«Ev»,*i\CJ)  o  0 1  0  1 Q  o 
kiRITE  16 1 1 )  T  » ( X(R,  I>  ,1=1, NE1 
1  FORMAT  <bX, 10H0.6) 

1.0.0  CuNTiNtLE 

RETURN 

END 
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Suokuu  1  1  NE  u  A  C  3  S  l  N ) 
i  i\  T  t  u  t  K  I  i  l  1 C  i  j  I  )  J  i  K.  t  N  »  L  » ivi  t  H  » M  M  »  1 
ricAL^b  A  (  1 0  t  10  j  >  rs  C  i  0  i  »  X  l  1  0  J  t  S  >  U 
HlAL  X  X  l 10) 

L L N /Ab//  A  ,k/ ABE/ .X A 

C 

C 

C  The  PURPOSE  Ur  IhlS  SUdKGUT  INb  lb  To  SOlVE  THE 

C  SYSTEM  OF  LINEAR  ALoEbkAlC  buUAl  IONS  ASSUCIAT  EU 

C  k  I  T  H  THE  LINEAR  LEAST  3  U  U  A  R  b  S  PkUbLEM.  THE  GAUSS 

L_ELIMINATION  PkUCLLDURE  IS  USED  y»ITH  THE  PIVuT 

L  ELEMENT  BEING  THE  LARGEST  IN  THE  UNREDUCED  PART 

C  LF  THE  MATRIX. 

L 

M=N-  I 

no  to  i  - 1 .  n 

10  II ( I )=I 

L)U  Ii  J=  1 »  M 

D L  12  I=J  #N 
DG  12  K  =  J  1 1\ 

_ n=oABS  t  A  l  I  »  J  i  I _ 

IF(G»lE»S)  GOTO  12 
S  =  U 

T  =  N 

lb  CONTINUE 

IFtL.EC.J)  GCTG  19 

DG  1 A  I  =  j  » N 
S=A(L,  I  ) 

ail,o=a u ,  n 

1 A  A  I J  »  i  )  =  b 
3  =  k ( L  > 

_ K  (  L  )  =  R  (  J  ) _ 

K  (  J  )  =  S 

I  ■)  I F ( K  «  EC . J i  Gu  TO  13 
DU  2J3  i  =  J,N 
S  =  A {  I  , T ) 

A(l»T)-A(i»J) 

_ Z‘j  A  (  1  t  J  )  =  S _ 

I  =  1  I  (  T  ) 

I I  t  T  I  =  1 1  <  J  ) 

1 1  l  J  )  =  i 

13  IFCGABSI AIJ*J) i.GT.l. 00-35)  ouTG  13 
IFIJ.EC.l  )  GLTu  13 
_ wKllc(fc»3) _ 

STOP 

13  M  M  =  J  + 1 
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1  1 

lit  t 

)  * 

l  t  / 

I E  (  L  A  b  $  (  A  {  I »  J  /  )  .LT.l.QD-33)  buTU  11 

- S~A.(.Jt  >1  i  /.All  ,  J  J 

A  {  I  ,  j ) =0 . 0 

CG  lc  K=IV^,K 

^G_AU^lUjiAXJ_LiYj-S»A  (  1  .Ki 

M  I  J  -  K  i  j  )  -  S  *  k  {  i  i 
i  1  CUi'J  r  I N  UE 
u^i  LI  in  —  i|i'( 

I  =  N  +  l  -  K 

S=0 . 0 

- LEI  I  .Fii.lSi  >  ftr-m  17 _ 

F  M  =  i  + 1 

DC  18  ,  is 

lo  S=S ±  A  (  1  .«  J 1  *  X I  .J  l 
17  XI  i )=( Ki  I  ) -b  J /A(  I  ,  n 

0U  2  i  i  =  1  ,  t\ 

_ 1=1111) _ _ _ _ 

IE(  i  .  L  G  •  K  j  bCFG  21 
S=X { K ) 

a  1 1\  J  =  XJL..L) _ „ _ _ _ 

X  {  1  }  =  8 

in  i  )  =  ii  (K) 

- U  (K  )  =  * _ _ _ 

21  CUMINUfc 
GO  23  1  =  1  ,N 
23  X  X 11 1= x  i  ii 

3  FflRWAT { 16HJMATRIX  SINGULAR) 

KG  TURN 
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oUdKLUTIM-  LEAST  ( l\  ,  M  ) 

!.N.ILG  ££.  I  *  j  » i\ ,  >  ,  ;> 

KtAL^o  A( 500, 101 ,B( 10,101 ,RC 500) ,CC 10) fSUMA, SUMS 

REAL  X(IO) 

- C  L F  F  u  N _ /  Ac2/t,C/Ab6/X/Atjl/A»K _ _ 

c 

C 

t  Tilt.  HURHCSt  C £  IHJLS  SUoKOuTiNt  LS  Tu  nbEbCt  THE 

C  l c m  j I  SQUARE  PROBLEM  7 J  A  SYSTEM  OF  LINEAR 

C  ALbtfcKAiC  EQUATIONS. 

L 


I F ( M  . G  T  .  N  J  uu  TO  1 A 
IFiM.EQ.MGGTG  17 
w  R  i  )  t  (  c  »  2  J 

2  FC  R  F  AT (18HJ INSUFFICIENT  DATA) 

ST  OF _ 

17  DC  1 5  I  =  1 »  N 
Cli)=R(IJ 
DC  15  J^1,N 
Id  l(1,o)-a!1,J) 

GCTC  16 

14  DO  10  I =1  .  N _ 

00  12  J  —  I  j  i\ 

S  U  F  A = G  .  0 
DC  U  0^1, 

i  1  SUMAwSUMA+A  IS»I1*AIS, J) 

D(  1  »  j  )  =  S  Uv  A 

Id  b  (  j  t  I  )  =  S  C  Fi  A _ 

SUFA=0 .0 
CO  Id  5  =  1 »  H 

ID  S  0  F  A  -  S  0  Fi  A  A  {  S  »  L  J  ^  K  IS) 

10  C  C  i  )  =  S  UM  A 
16  CALL  GAUSS(N) 

_ RE  I  L k N _ 

END 
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C 

C 

JL. 

C 

0 

_L_ 

C 

c 

JL. 

C 

l 

c 

c 

c 

X__. 

c 

c 


SUBKGU1  I  N  E  SLCFt£CR,T  ,  I  JACK) 

fi.£,4L  10. )  >  o & 1  LQiIQ I  iG i  10 )  iT,C(  1 1 ,  l  j  }  ,XQ( 3 

C  ,  1 C  ) 

INTEGER  f<  ,NE  ,NK,  IMG  »  I  I ,  NOS  »  1  J  ACK 

Lbh^GN  /  A  Li  3/  L  ,  NE  t  NK  t  NO /Abb/  JA  ,  jt\  ,  X  .b/  Auh/XC*  1  1  ,  NtJ.S _ 


1  h  1  ^ . S.U.URCUT1N.E . HU  ST  BE  SUPPLIED . ..ALCNG  WITH  THE _ 

FOR  THE  PARTICULAR  PROBLEM*  THE  HjHPu^L  OF  THIS 
SUBROUTINE  IS  Tu  EVALUATE  ThE  FUNCTIONS  AND  THE 

JACUbiAN  f^ATRj  C  1  E  S  . _ 3  HE  FORMAT  lF  TtiiS  MipKuUT  INfc 

MUST  bt  AS  SHOWN  IN  THIS  EXAMPLE. 


lhL  SUBROUTINE*  REAL*  INTEGER  AND  COMMON  STATEMENTS 
FLS1  bE  SOPPLiEu  EXACTLY  AS  AbGVE. 


AUDIT  1CNAL  KLALt  INTECER  UK  uTHER  TYPE  STATEMENTS 
IvmY  hi.  SUPPLIED  IMMEDIATELY  BELOW  I F  NECESS  ARYt 
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c 
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c_ 
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-EX A  1 — X.l,X2,MtK2,Kb  ,KEi,l\E2«Rr-T 

IM  EGER  I  i  j 


uflCK  Tu  DETERMINE  IF  3 hl  STATE  VARIABLE  BOUND  IS 
LXlLtuCb.  IF  SO  THE  EXECUTION  OF  THE  PROGRAM  IS 

IXjliiJ-jVA  X  l-  . _ 3hi_KE  ARE  SuVLrsAU  kEa^UnS  Fuk  FAILURE 

AT  THIS  flint,  Sure  OF  WHICH  ARE: 

LJJHJE  DIFFERENT  I AL  EQUAT  IONS  ARE  NOT  PROPERLY 
NORMALIZED. 

_ l  NUMEkIlAL  INTEGRATION  PROLb  EUURE  IS 

uNslAbLt  FUR  IhlS  SYSTEM  wiTh  TnE  olVeN  NO. 

OF  INTEukATIuN  POINTS  N  C  (  I  ,  J  )  . 


3  CHMAX  IS  TOO  LARGE  AND  THEREFORE,  THE  STEP 
bluE  RMAXI  Is  AL  SL  TOO  LARGE. 


DU  It  1*1 ,  K  E 

IF (AdStX IF, 1 j) .Lt. 10.0 )  GOTO  16 

WR I T  E 16  *  1 1  i 

STOP 

lfc  CONTINUE _ 

17  FORMA  1  (• OSTA It  VARIAbLE  BCUNL  EXCEEDED') 
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APPENDIX  B 

The  1-Butene  Problem 

For  Example  Two,  exact  data  was  generated  from 
the  model 


dx. 


dt 


=  -a, (x,  -  k  , x„ )  -a 

i  ±  el  2 


al  -ke2(1'xl-x2 


0 


dx2 

dt 


al  (xl  kelX2}  a3  \^2  ke3(1  X1  X2^ 

(B-l) 


A  normally  distributed  error  was  then  obtained  from  a 
random  number  generator.  Then  each  exact  data  point  was 
perturbed  by  setting 


x 


ijk 


=  x 


ijk 


(1  + 


0£  .  ) 
i  jk 


i  =  1 


=  1 


k  =  1 


with  representing  an  experimental  observation  of 

x.  which  is  the  actual  value  of  the  ith  state  variable 
13k 

for  the  jth  experiment  at  the  kth  discrete  value  of  time. 
The  relative  error  ae^^.  f°r  this  point  has  a  standard 
deivation  of  a.  Three  values  of  a  (0.1,  0.01,  0.001)  were 
used  to  perturbate  the  data  generated  from  equation  B-l 


B-2 


with  the  parameters  and  equilibrium  constants  listed 
below 


al  = 

10.344 

it 

i — 1 
0) 

0.4469 

a2 

3.724 

k  0 
e2 

0.2685 

(B-3 ) 

a3  = 

5.616 

k  o  = 
e3 

0.6002 

and  the  initial  conditions  as  given  in  Table  B-l. 


TABLE  B-l 

Initial  Conditions  for  Example  Two 


Experiment  Number 

(0) 

x2(0) 

1 

1.0000 

0.0000 

2 

0.0000 

1.0000 

3 

0.0000 

0.0000 

4 

0.2400 

0.7600 

5 

0.3492 

0.6508 

6 

0.0000 

0.4937 

7 

0.4130 

0.0000 

Several  different  attempts  were  made  to  determine 
the  parameter  vector  from  the  data.  Different  combina¬ 
tions  of  experiments  were  used  to  test  the  sensitivity 
of  the  method  to  the  choice  of  experiment  as  well  as  to 
experimental  error.  The  results  of  these  attempts  are 
presented  in  the  following  tables.  Each  attempt  was  made 


• 

B-3 


with  the  initial  guess  of 


a 


(0) 


1.0 

1.0 


1.0 


(B-4 ) 


for  the  parameter  vector.  This  guess  was  sufficiently 
close  to  the  answer  to  use  quasilinearization  directly. 

Tables  B-2,  B-3,  and  B-4  contain  the  convergent 
sequences  of  parameters  when  all  the  seven  experiments 
are  included. 


TABLE  B-2 
1-Butene  Problem 
All  Seven  Data  Sets  a  =  0.1 


Iteration  a^  a^  a^ 


1.000 

1.000 

1 .000 

2.389 

2.417 

2.523 

5.820 

3.727 

4.449 

9.937 

3.790 

5.380 

10.845 

3.662 

5.512 

10 . 943 

3.616 

5.540 

10 .929 

3.619 

5.540 

10.930 

3.619 

5 . 540 

7 


B-4 


TABLE  B-3 
1-Butene  Problem 
All  Seven  Data  Sets  a  =  0.01 


Iteration 

al 

a2 

a3 

0 

1.000 

1.000 

1.000 

1 

2.383 

2.448 

2 . 546 

2 

5.584 

3 .885 

4.532 

3 

0.409 

3.885 

5.482 

4 

10.381 

3.721 

5.605 

5 

10.404 

3.713 

5.609 

6 

10.404 

3.713 

5.609 

All 

Iteration 

TABLE 

1-Butene 

Seven  Data 

al 

B-4 

Problem 

Sets  a  =  0.001 

a2 

a3 

1.000 

1.000 

1 . 000 

2.383 

2.451 

2 . 549 

t .  t6  0 

3.901 

4.541 

0.357 

3.894 

5.491 

10.331 

3.728 

5.613 

10.350 

3.722 

5.615 

5 


■ 


B-5 


In  the  following  three  tables,  the  results  of 
using  the  data  given  in  the  two  eigenvector  directions 
are  presented.  Experiments  five  and  six  are  in  the  two 
respective  eigenvector  directions. 


TABLE  B-5 

1-Butene  Problem 
Data  Sets  Five  and  Six  a  =  0.1 


Iteration 


al 

a2 

a3 

1.000 

1.000 

1.000 

2.608 

2.407 

2.515 

6.801 

3.515 

4.424 

10.751 

2.948 

5.615 

12.464 

2.293 

6.088 

12.111 

2 . 225 

6.135 

12.247 

2.242 

6.146 

12 . 222 

2.242 

6.146 

12.223 

2 . 242 

6.146 

12.223 

2 . 242 

6.146 

0 

1 

2 

3 

4 

5 

6 

7 

8 


9 


' 


B-6 


TABLE  B-6 

1-Butene  Problem 
Data  Sets  Five  and  Six  a  =  0.01 


Iterations 

ai 

a2 

a3 

0 

1.000 

1.000 

1.000 

1 

2.407 

2.434 

2.547 

2 

5.682 

3.844 

4.537 

3 

9.540 

3.789 

5 . 523 

4 

10.504 

3.598 

5.653 

5 

10.513 

3.588 

5.658 

6 

10.513 

3.588 

5.658 

TABLE 

1-Butene 

Data  Sets  Five 

Iteration  a^ 

B-7 

Problem 

and  Six  a  = 

a2 

0.001 

a3 

0 

1.000 

1.000 

1.000 

1 

2.387 

2.436 

2 . 550 

2 

5.567 

3 . 877 

4.548 

3 

0.362 

3.877 

5.510 

4 

10.343 

3.714 

5.618 

5 

10.361 

3.710 

5.620 

Using  the  data  from  the  first  experiment,  the  results  in 
Tables  B-8,  B-9,  and  B-10  were  obtained. 


• 

B-7 


TABLE  B-8 
1-Butene  Problem 
Data  Set  One  a  =  0.1 


Iteration 

ai 

a2 

a3 

0 

1.000 

1.000 

1.000 

1 

4.670 

0.198 

34.811 

2 

12.749 

-1.617 

23.507 

3 

5.855 

7.626 

-2.158 

4 

9.454 

5.182 

3.252 

5 

8.796 

5.725 

1.782 

6 

9.744 

4.806 

3.515 

7 

9.361 

5.173 

2.815 

8 

9.673 

4.887 

3.337 

9 

9.519 

5.028 

3.074 

10 

9.615 

4.941 

3.235 

complete  convergence  not  obtained 


B-8 


TABLE  B-9 
1-Butene  Problem 


Data  Set  One  a  =  0.01 


Iteration 

ai 

a2 

a3 

0 

1.000 

1.000 

1.000 

1 

4 .321 

5.667 

30.005 

2 

12.267 

-1.394 

24.771 

3 

0 . 281 

4.514 

2 . 203 

4 

10.275 

3.873 

5 . 511 

5 

10.177 

3.923 

5 . 203 

6 

10.167 

3.932 

5.188 

7 

10.167 

3.932 

4.189 

TABLE 

1-Butene 

Data  Set  One 

Iteration  a^ 

B- 10 

Problem 

a  =  0.001 

a2 

a3 

0 

1.000 

1.000 

1.000 

1 

4.287 

0.604 

29.525 

2 

12.224 

-1.376 

24.913 

3 

9.784 

4.061 

3.017 

4 

10.453 

3.659 

5 . 894 

5 

10 . 325 

3.745 

5.570 

6 

10.324 

3.746 

5.569 

For  data  from  the  second  experiment  the  results 


in  Tables  B-ll,  B-12,  and  B-13  wepe  obtained. 


B-9 


TABLE  B-ll 

1-Butene  Problem 
Data  Set  Two  a  =  0.01 


Iteration 


0 

1.000 

1.000 

1.000 

1 

4.403 

31.596 

1.686 

2 

15.161 

31.968 

1.825 

3 

7.605 

8.436 

6.205 

4 

11.859 

10.018 

5 .230 

5 

14.661 

14.493 

4.516 

6 

18.369 

21.451 

3 . 464 

7 

28.487 

38.500 

0.678 

8 

775.187 

1092.55 

193.895 

The  integration 

of  fundamental  matrix  became 

unstab 

a  solution  was  : 

not  obtained 

• 

TABLE 

B-12 

1-Butene 

Problem 

Data  Set  Two  a  =  0.001 

Iteration 

ai 

a2 

a3 

0 

1.000 

1.000 

1 .000 

1 

4.614 

32.859 

1.622 

2 

11.088 

23.740 

3.063 

3 

0.733 

-0 . 267 

5.722 

4 

10.404 

4.016 

5 . 600 

5 

10 . 355 

3.758 

5.614 

6 

10.355 

3.758 

5.614 

B-10 


Attempts  were  also  made  to  identify  the  parameters 
from  experiment  three,  from  experiment  five,  and  from 
experiment  six.  All  three  of  these  attempts  were  com¬ 
plete  failures.  The  results  diverged  and  were  out  of 
bound  within  a  few  iterations. 

To  test  the  sensitivity  of  the  method  to  error  in 
the  initial  conditions,  perturbations  were  introduced 
into  the  initial  conditions  for  experiments  one,  two  and 
three.  Attempts  were  then  made  to  determine  the  para¬ 
meters  using  the  first  three  data  sets  with  standard 
deviations  of  0.001  and  0.1. 

First  with  the  exact  initial  conditions  of 


a1  (0) 


1.0 

0.0 


X2(0) 


0.0 

1.0 


x3(0) 


0 . 0 
0.0 


(B-5) 


the  solutions 


a 


10.346 

3.727 

5.631 


and 


a 


10.638 

3.930 

5.471 


(B-6 ) 


(B-7 ) 


were  found  for  standard  deviations  of  0.001  and  0.1 
respectively . 


B-ll 


Perturbating  the  initial  condition  to 


0.99 

0.01 

0.01 

II 

o 

1 — 1 
x\ 

1 

i — 1 

O 

• 

o 

7  X2  (°)  = 

_  0 . 9  9  _ 

;  x3 (0)  = 

1 

1 — 1 
O 

« 

O 

_ 1 

produced  the  results 


a 


and 


a 


10.177 

3.724 

5.545 

10.479 

3.921 

5.389 


(B-8 ) 


(B-9 ) 


(B-10 ) 


respectively  for  the  two  standard  deviations  of  0.001 
and  0.1. 

With  the  initial  conditions  perturbated  to 


xx  (O) 


0.9 

1.0 


x2(0) 


0.1 

0.9 


—3  (0) 


0.1 

0.1 


(B-ll) 


the  results  for  the  standard  deviation  of  Q.001  were 


8.487 


a  •= 


3.968 


4.876 


(B-12 ) 


and  for  the  standard  deviation  of  0.1  they  were 


' 

B-12 


a 


8.866 

4.143 

4.759 


(B-13 ) 


The  input  data  which  was  used  to  produce  the 
results  in  Tables  B-2,  B-3,  and  B-4  is  presented  immedi¬ 
ately  following  the  Fortran  IV  listing  of  the  SLOPE 
subroutine  required  for  this  problem.  The  input  data  used 
for  the  other  results  can  be  formed  by  deleting  the  data 
sets  for  the  appropriate  experiments. 


B-l  3 


SU8RCUT I NE  S LCPE ( K  ,  T  ,  I JACK ) 

6  L  U  P  E  GUBKCuTINE  FUR  frtc  1  —  oUTENE  PROBLEM 

RERL  X( l3,10 ) , JX< 10,10) » JK( 10, 10) ,G( 10) ,T,C(llt 10)  ,  XuC  3 
1C, 10) 

_ INTEGER  R  ,NE  .Ni\,NC,  I  1  ,NDS  .  I  JACK _ 

COMMON  /AB3/C,N£,NK,N0/AB3/JX,JKfX,G/AB8/Xu,  i  I  ,  NDS 
REAL  XI »X2»K1,K2,K3,KE1»KE2»NE3 
INTECER  1,J 
CO  16  1*1, NE 

IF i AB S ( X { R ,  1  )  )  .LE.  10  .0  )  G  0  1  0  16 
ftRlTElo,  17) _ _ 

STOP 

16  CONTINUE 

17  FUkMAT { ' OS  1A  IE  VARIABLE  SOUND  EXCEEDED1) 

X  1—  X  (  R  ,  1  ) 

X2  =  X {  R  ,  2  ) 

_ K.1  =  C<  R  ,  1  ) _ _ 

K2=C l R , 2  ) 

K3=C{ R,3) 
mim n,i2) , Xu 

11  DC  13  1  =  1,  NE 
G ( I  )  -0  .  G 

DO  14  J  =  1  .Nr _  _ 

14  Ja( 1 ,  J  )  =G  .  0 
DO  13  J= 1 , NK 
13  JKU  ,  J)  =  C.O 
Nu=2 

KE1=0.4469 

RL2= 0.2683 _ 

KE  3=  C . 6CC2 

12  G< l)=-Kl^(Xl-KEl^X2)-K2*(Xl-KE2^I 1.0-X2-X1) ) 

G{ 2)=K1*(X1-KE1*X2)-K3*i  X2~K£3*l 1 »Q-X2~X1)  ) 

IF  C I JACK.LT* 0)  GOTO  13 

J  X ( 1 , 1 ) =-Kl-K2-KE2*K2 

JX(l,2)  =  Rl*KEl-K2»KE2 _ 

JX( 2 , 1 )=K1-K3*KE3 

J X  (  2,2)  =-Kl*KEl-K3-K3*KE3 

J  K.  I  1  ,  li-~Xi  +  KEl>!tX2 

JR{ 1  ,2  )=-Xl  +  KE25M  1.0-X2-A1) 

JK(  2, 1 )  =  X 1-K E  1  *X2 

JK  (  2 ,  3 j  =  -X2  +  KE3*(  1  « C-X2-X 1 J _ 

13  RETURN 
END 
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APPENDIX  C 

Non-Isothermal  Pyrolysis  of  Propane 

To  model  the  non-isothermal  pyrolysis  of  pro¬ 
pane  in  a  tubular  reactor,  the  following  equation  is 
sufficient. 


df  s 

-  E  -i 

A 

P  (L ) 

a 

1  -  f 

dL  „  exp 
FxlO 

R'T(L)_ 

_RT  (L)_ 

N 

L  +  p2  +  f 

(C-l) 

To  determine  the  parameters  A,  E,  and  a,  data 

from  sixteen  different  experiments  was  available. 

The  conversion  at  the  outlet  of  the  reactor  f (L  ) 

m 

was  measured  for  each  experiment.  The  temperature  was 

measured  at  twenty-seven  equally  spaced  points,  so  that 

intermediate  points  could  be  accurately  established  by 

linear  interpolation.  The  pressure  was  measured  at  both 

the  inlet  and  outlet  of  the  reactor.  Intermediate  values 

were  determined  by  assuming  that  pressure  was  a  linear 

function  of  L.  The  inlet  feed  rate  of  propane  F  and  of 

inlet  feed  rate  of  inerts  N  were  also  observed.  All  of 

o 

this  information  is  presented  in  Table  C-l. 

In  addition  to  this  information,  the  physical 
dimensions  of  the  reactor  must  be  specified.  The  cross 
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2 

sectional  area  s  is  0.0742242  cm  .  The  length  L  is 

J  m 

69.0118  cm. 

The  constants  R  and  R'  in  equation  (C-l)  are 
the  ideal  gas  constants  with  the  appropriate  units. 
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TABLE  C-l* 

Data  for  Propane  Pyrolysis 


Data  Set  No. 

1 

2 

3 

4 

5 

6 

7 

8 

FX104 (2rmoIe) 

8.80 

5.30 

4.26 

1.93 

1.85 

2.69 

3.57 

1.72 

N^04<2!f±e> 

41.5 

28.1 

17.5 

46.4 

31.6 

30.9 

29.9 

19.9 

Pi  (mmHg) 

974 

896 

834 

986 

915 

91b 

900 

839 

Pe (mmHg) 

742 

742 

742 

740 

744 

734 

734 

739 

Z  (L  ) 
m 

0.666  0 

0.5890 

0.6685 

0.6173 

0.6090 

0.6318 

0.6992 

0.6952 

Temp,  at 

340 

290 

280 

364 

462 

486 

443 

510 

Equally 

631 

545 

584 

646 

761 

792 

746 

842 

Spaced  Points 

928 

820 

875 

962 

1055 

1087 

1029 

1116 

(6F) 

1082 

1057 

1123 

1182 

1242 

1273 

1224 

1288 

1166 

1190 

1264 

1300 

1359 

1382 

1335 

1397 

1218 

1286 

1364 

1374 

1447 

1460 

1417 

1474 

1279 

1363 

1444 

1444 

1518 

1523 

1486 

1530 

1376 

1442 

1507 

1520 

1575 

1578 

1549 

1567 

1430 

1504 

1552 

1581 

1612 

1617 

1593 

1584 

1437 

1551 

1581 

1617 

1634 

1637 

1618 

1589 

1456 

1574 

1593 

1636 

1637 

1641 

1628 

1589 

1499 

1598 

1593 

1656 

1637 

1641 

1629 

1577 

1583 

1617 

1591 

1682 

1633 

1640 

1629 

1559 

1621 

1635 

1588 

1693 

1622 

1631 

1622 

1533 

1638 

1634 

1578 

1690 

1602 

1611 

1608 

1497 

1647 

1637 

1558 

1678 

1578 

1587 

1585 

1454 

1650 

1627 

1534 

1651 

1543 

1551 

1551 

1402 

1647 

1610 

1499 

1619 

1503 

1512 

1512 

1350 

1623 

1582 

1462 

1577 

1449 

1461 

1463 

1288 

1589 

1548 

1409 

1525 

1387 

1403 

1407 

1221 

1542 

1501 

1351 

1454 

1308 

1324 

1329 

1133 

1474 

1438 

1272 

1367 

1213 

1228 

1235 

1020 

1412 

1349 

1169 

1277 

1123 

1129 

1141 

920 

1238 

1262 

1082 

1152 

964 

964 

986 

725 

1070 

1033 

841 

939 

744 

747 

779 

517 

873 

810 

604 

738 

526 

543 

573 

336 

670 

595 

405 

560 

360 

390 

415 

205 

*  This  dat 

a  has 

been  t 

.aken  f 

'rom  Kershenbaum ' s 

Ph.D. 

dissertation  (University  of  Michigan  1964) 

: 

TABLE  C-l  (Continued) 
Data  for  Propane  Pyrolysis 


Data  Set  No. 

9 

10 

11 

12 

13 

14 

15 

16 

Fxl04(g-mOle) 

4  sec 

1.69 

2.33 

0.93 

0.80 

0.82 

1.01 

1.31 

1.42 

n  xicr 
o 

19.1 

27.4 

53.3 

55.3 

37.0 

50.7 

40.3 

25.2 

Pi  (mmHg ) 

839 

890 

1013 

994 

1014 

1001 

927 

854 

Pe (mmHg) 

739 

728 

727 

727 

740 

740 

740 

740 

Z  (L  ) 
m 

0.7909 

.4771 

.  6391 

.3740 

.3876 

.5592 

.7197 

.8022 

Temp .  at 

233 

321 

182 

342 

396 

522 

516 

545 

Equally 

543 

713 

449 

706 

718 

805 

792 

836 

Spaced  Points 

858 

1063 

798 

1057 

972 

1039 

1028 

1074 

(°F ) 

1117 

1316 

1107 

1300 

1117 

1151 

1154 

1217 

1273 

1459 

1268 

1431 

1188 

1217 

1234 

1312 

1365 

1550 

1364 

1504 

1226 

1292 

1302 

1384 

1433 

1618 

1433 

1569 

1408 

1509 

1473 

1476 

1479 

1671 

1507 

1642 

1609 

1619 

1579 

1537 

1506 

1710 

1588 

1720 

1684 

1677 

1629 

1558 

1514 

1720 

1637 

1749 

1726 

1708 

1654 

1558 

1514 

1719 

1663 

1767 

1751 

1727 

1668 

1553 

1511 

1716 

1687 

1798 

1772 

1741 

1669 

1543 

1499 

1710 

1711 

1832 

1789 

1751 

1667 

1526 

1477 

1696 

1739 

1843 

1797 

1751 

1659 

1502 

1449 

1683 

1741 

1840 

1797 

1751 

1642 

1472 

1411 

1654 

1740 

1824 

1793 

1739 

1620 

1437 

1373 

1624 

1723 

1801 

1783 

1719 

1592 

1397 

1325 

1583 

1692 

1771 

1766 

1691 

1554 

]  349 

1276 

1541 

1656 

1728 

1737 

1696 

1507 

1297 

1216 

1479 

1602 

1677 

1697 

1609 

1454 

1235 

1149 

1406 

1539 

1605 

1638 

1542 

1381 

1159 

1063 

1312 

1461 

1517 

1564 

1468 

1303 

1073 

953 

1187 

1358 

1415 

1483 

1363 

1207 

971 

862 

1081 

1281 

1296 

1273 

1101 

963 

742 

676 

829 

1073 

1050 

929 

786 

642 

454 

489 

596 

879 

845 

682 

583 

455 

293 

323 

405 

686 

655 

530 

420 

334 

205 

• 

' 
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Two  problems  were  considered.  In  the  first,  a 
was  fixed  at  one  and  the  other  two  parameters  were 
determined  such  that  the  least  square  error  was  mini¬ 
mized.  Secondly,  all  three  parameters  were  determined 
simultaneously . 

The  maximum  step  size  S  ,  the  maximum  allowable 

^  max 

relative  change  in  the  parameters  for  one  iteration  cmaxr 

and  the  minimum  allowable  relative  change  in  the  parameters 

for  one  step  C  .  which  were  used  with  the  different  initial 
^  min 

guesses  for  A  and  E/R  are  specified  in  Table  C-2.  The  360/67 
machine  time  is  also  given  in  this  table. 


TABLE  C-2 

Control  Data  for  First  Propane  Problem 


Initial 

A 

Guess 

E/R 

S 

max 

C 

max 

C  . 
mm 

IBM  360/67 

Time  (min) 

35.0 

26.0 

2xl0"3 

2.0 

0 . 5 

3.90 

10.0 

10.0 

-4 

1x10  q 

5.0 

0.5 

4.56 

30.0 

30.0 

-4 

1x10  q 

5.0 

0.5 

4.34 

40.0 

40.0 

-4 

1x10  4 

5.0 

0.5 

4.44 

0.0 

0.0 

-4 

1x10 

5.0 

0.5 

7.73 

50.0 

50.0 

lx.10"4 

5.0 

0.5 

4.03 

35.4 

26.22 

lxlO-4 

5.0 

0.5 

3.98 

The  convergent  sequences  of  parameters  which  evolved 
from  each  of  the  initial  guesses  in  Table  C-2  are  shown  in 
Tables  C-3,  C-4,  C-5,  C-6,  C-7 ,  C-8,  and  C-9 . 

TABLE  C-3 


First  Propane  Problem  First  Initial  Guess 


Step 

No. 

Sum  of 

A 

E/R 

Errors 

Sauared 

Comment 

0 

35.00 

26.00 

.  4456 

Begin  Phase  I 

1 

34.92 

25.90 

.4390 

2 

34.77 

25.71 

.4260 

3 

34.81 

25.33 

.4005 

4 

33.94 

24.63 

.  3525 

5 

32.99 

23.41 

.  2679 

6 

31.54 

21.52 

.1424 

7 

29.80 

19.20 

.0382 

Begin  Phase  II 

8 

29.46 

18.73 

.0335 

Solution 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 
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TABLE  C-4 

First  Propane  Problem  Second  Initial  Guess 


A 

E/R 

Sum  of 
Errors 
Squared 

Comment 

10.00 

10.00 

2.535 

Begin  Phase  I 

15.41 

14.98 

2.535 

16.24 

15.40 

2.535 

17.08 

15.73 

2.535 

17.88 

15.94 

2.534 

18.63 

16.06 

2.534 

19.37 

16.13 

2.532 

20.08 

16.17 

2.530 

20.79 

16.19 

2.525 

21.49 

16.20 

2.515 

22.19 

16 . 21 

2.494 

22.90 

16 . 23 

2.454 

23.62 

16.25 

2.375 

24.36 

16.31 

2 . 219 

25.17 

16.41 

1.925 

26 . 10 

16.65 

1.399 

27.33 

17.20 

0.6114 

Begin  Phase  II 

29.45 

18.72 

0.0345 

Solution 

IN  w  • 
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16 
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TABLE  C-5 

First  Propane  Problem  Third  Initial  Guess 


A 

E/R 

Sum  of 
Errors 
Squared 

Comment 

30.00 

30.00 

2.534 

Begin  Phase 

2914 

28.86 

2.534 

27.79 

26.06 

2.534 

26.02 

24.66 

2.533 

24.19 

22.06 

2.533 

22.78 

10.83 

2.532 

22.00 

18.26 

2.529 

21.81 

17.30 

2.524 

22.02 

16.78 

2.514 

22.46 

16.50 

2.494 

23.03 

16 . 37 

2.453 

23.68 

16.32 

2.374 

24.40 

16 . 34 

2.219 

25.18 

16.43 

1.924 

26 . 10 

16.65 

1.398 

27.33 

17 . 20 

0.6109 

Begin  Phase 

29.45 

18.72 

0.0336 

Solution 

LN  'w'  • 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 


C-9 


TABLE  C-6 


First  Propane  Problem  Fourth  Initial  Guess 


A 

E/R 

Sum  of 
Errors 
Squared 

40.00 

40.00 

2.530 

39.73 

39.65 

2.530 

39.20 

38.98 

2.529 

38.23 

37.73 

2.529 

36.56 

35.58 

2.529 

34.03 

32.30 

2.527 

30.88 

28.16 

2.525 

27.86 

24.07 

2.520 

25.67 

20.86 

2.510 

24.48 

18.76 

2.489 

24.10 

17.56 

2.449 

24.22 

16.93 

2.370 

24.66 

16.64 

2.214 

25.31 

16.57 

1.920 

26.16 

16.72 

1.394 

27.36 

17.23 

0.6081 

29.45 

18.72 

0.0336 

Comment 


Begin  Phase  I 


Begin  Phase  II 


Solution 


. 

Step 

No. 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 
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TABLE  C-7 


First  Propane  Problem  Fifth  Initial  Guess 


Sum  of 

A 

E/R 

Errors 

Scruared 

0.00 

0.00 

2 , 535 

0.00 

0.00 

2.535 

0.00 

0.00 

2.535 

0.00 

0.00 

2.535 

0.00 

0.00 

2 . 535 

0.00 

0.00 

2.535 

0.00 

0.00 

2.535 

0 . 00 

0.00 

2.535 

0 . 26 

0.19 

2 . 535 

0 . 80 

0.62 

2.535 

1.48 

1.20 

2 . 535 

2.06 

1.71 

2.535 

2.67 

2.27 

2.535 

3.20 

2.76 

2.535 

4.08 

3.60 

2.535 

4 . 86 

4.36 

2.535 

6 . 08 

5.55 

2.535 

7.00 

6.46 

2.535 

8.48 

7.91 

2.535 

10.43 

9.79 

2.535 

12.56 

11.78 

2.535 

14.52 

13.46 

2.535 

16.11 

14.63 

2 .535 

17.36 

15.35 

2.534 

18 . 36 

15.75 

2.534 

19.23 

15.97 

2.532 

20.01 

16.09 

2.530 

20.75 

16.15 

2 . 525 

21.47 

16.18 

2.515 

22.18 

16 . 20 

2.494 

22.89 

16.22 

2.454 

23.62 

16 . 25 

2.375 

24 . 36 

16 . 30 

2.219 

25.17 

16.41 

1.925 

26.10 

16.65 

1.399 

27.33 

17 .20 

0.6114 

29.45 

18.72 

0.0336 

Comment 


Begin  Phase  I 

s-  -*  s  /2 

max  max,- 

s  -*  s  /2 

max  maxi, 

s  -*  s  /2 

max  maxi, 

s  s  /2 

max  max^ 

s  ■*  s  /2 

max  max,, 

s  -*  s  /2 

max  maxi, 

s  +  s  /2 

max  max' 


Begin  Phase  II 
Solution 


■ 


0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 
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TABLE  C-8 

First  Propane  Problem  Sixth  Initial  Guess 


A 

E/R 

Sum  of 
Errors 
Squared 

Comment 

50.00 

50.00 

2.502 

Begin  Phase  I 

49.94 

49.92 

2.501 

49.82 

49.77 

2.501 

49.59 

49.48 

2 . 501 

49.14 

48.90 

2.501 

48 . 27 

47.79 

2.499 

46.66 

45.74 

2.497 

43.92 

42.24 

2.492 

39.88 

37 . 06 

2.483 

35.15 

30.93 

2.463 

30.97 

25.37 

2.433 

28.17 

21.39 

2 . 343 

26.75 

18.99 

2.187 

26.35 

17.73 

1.892 

26 . 64 

17.24 

1.367 

27.54 

17.43 

0.488 

Begin  Phase  II 

29.45 

18.72 

0.0336 

Solution 

N— '  • 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 
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TABLE  C-9 

First  Propane  Problem  Seventh  Initial  Guess 


A 

E/R 

Sum  of 
Errors 
Squared 

Comment 

35.40 

26.22 

0 . 2728 

Begin  Phase 

35.40 

26 .22 

0.2728 

35.40 

26.22 

0 . 2728 

35.40 

26.22 

0 . 2727 

35.39 

26.21 

0.2726 

35.59 

26 . 21 

0.2722 

35 . 38 

26.19 

0 . 2715 

35.34 

26.16 

0 . 2701 

35.02 

26.10 

0 . 2673 

35.20 

26.97 

0 . 2617 

35.01 

25.73 

0 . 2507 

34.63 

25.26 

0.2297 

33.93 

24 . 38 

0.1913 

32.71 

22.85 

0.1289 

30 . 89 

20.54 

0.0564 

Begin  Phase 

29.45 

18.72 

0 .0336 

Solution 

C-13 


For  the  second  case,  with  a  being  unknown,  the 
computer  time  required  was  considerably  larger,  as  shown 
in  Table  C-10.  Also  given  in  Table  C-10  are  the  three 
initial  guesses  which  were  used  and  the  corresponsing 
initial  control  variables. 

TABLE  C-10 


Control 

Data 

for  Second 

Propane 

Problem 

A 

Initial  Guess 

E/R  a 

S 

max 

C 

max 

C  . 
min 

IBM  360/67 
Time  (Min) 

18 

15 

1.0 

2xl0~3 

2.0 

0.5 

6.94 

18 

15 

0.5 

lxlO"4 

2 . 0 

0 . 5 

8 . 35 

18 

15 

1.5 

lxlO"7 

2.0 

0.5 

>  15 

The  convergent  sequence  of  parameters  for  each  of  the 
initial  guesses  given  in  Table  C-10  are  listed  in  Tables  C-ll, 
C-12,  and  C-13.  The  third  initial  guess  was  not  carried  all 
the  way  to  solution  because  of  the  large  amount  of  computing 


time  required. 


. 


Step 

No. 

1 

2 

3 

4 

5 

6 

7 

8 

9 


C-14 


TABLE  C-ll 


Second  Propane  Problem  First  Initial  Guess 


A 

E/R 

18.00 

15. 00 

25.87 

20.14 

26.00 

20.06 

26  17 

19.97 

26.49 

19.82 

27 . 12 

.19.62 

28.34 

19.55 

30.90 

20.24 

33.43 

21.48 

a 

Sum  of 
Errors 
Squared 

1.000 

2.533 

0.749 

0.9034 

0.7599 

0 . 8717 

0.7740 

0.8304 

0.800 

0.7513 

0.844 

0.6068 

0.916 

0.3711 

1.025 

0.944 

1.109 

0.0303 

Comment 


Begin  Phase  I 


Begin  Phase  II 


Solution 


IN  W  • 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 


C-15 


TABLE  C-12 


Second  Propane  Problem  Second  Initial  Guess 


A 

E/R 

a 

Sum  of 
Errors 
Squared 

Comment 

18.00 

15.00 

0.500 

0.7150 

Begin  Phase 

18.01 

15.00 

0.500 

0.7146 

18.02 

15.01 

0.501 

0.7138 

18.05 

15.02 

0.502 

0.7122 

18.12 

15.04 

0.505 

0.7090 

18.24 

15.89 

0.510 

0.7025 

18.47 

15.18 

0.520 

0.6897 

18.93 

15.34 

0.540 

0.6644 

10.78 

15.66 

0.576 

0.6155 

21.30 

16.20 

0.641 

0 . 5241 

23.80 

17.12 

0.748 

0.3671 

27.76 

18.70 

0.908 

0.1496 

Begin  Phase 

33.43 

21.48 

1.109 

0.0303 

Solution 

C-16 


TABLE  C-13 


Second  Propane  Problem  Third  Initial  Guess 


Step 

No. 

A 

E/R 

a 

Sum  of 
Errors 
Squared 

Comment 

0 

18.00 

15.00 

1.500 

2.535 

Begin  Phase 

1 

12.49 

13.12 

1.180 

2.535 

2 

11.22 

12.37 

1.111 

2.535 

3 

10.93 

12.68 

1.042 

2.535 

4 

11.91 

13.59 

1.022 

2.535 

5 

11.95 

13.61 

0.987 

2.535 

6 

12.69 

14.05 

0.972 

2.535 

7 

13.80 

14.64 

0.969 

2.535 

8 

14.84 

15.07 

0.971 

2 . 535 

9 

15.66 

15.26 

0.970 

2.535 

10 

16.45 

15.38 

0.970 

2.535 

11 

17.19 

15.45 

0.970 

2.534 

12 

17.92 

15.49 

0.970 

2 . 533 

13 

18.63 

15.52 

0.970 

2 .532 

14 

19.34 

15.53 

0.970 

2 . 528 

15 

20.04 

15.54 

0.970 

2.522 

16 

20.75 

15.55 

0.970 

2.509 

17 

21.47 

15.57 

0.971 

2.483 

18 

22.21 

15.61 

0.972 

2.432 

19 

22.99 

15.68 

0.973 

2 . 331 

20 

23.87 

15.82 

0.977 

2.135 

Stopped  aft' 

21 

24.45 

16.13 

0.984 

1.769 

22 

26.52 

16.82 

1.000 

1.141 

15  minutes 

The  required  SLOPE  subroutines 

for  the 

two  problems 

are 

listed  on  the 

following 

pages . 

017 


SUBROUTINE  oLCPE I R , T , I JALK ) 

U  SJL£L££  SUBROUTINE  FOR  THE  ERJ2M&E  EBSl&LEh  teUh  The 

C  CRQER  jH  REACTION  FIXED  AT  UNITY, 

REAL  X(10,10),JXtl0,10)iJK{10,10),G(10),T,CUl,10i  ,XOI  3 

_ 10*10) _ 

I N T t G E R  R,N£,NK»ISiA*Il,NUS,IJACK 

COMMON  /AB3/C,N£, NK» NA/AB5/JX, JK ,X,G/A68/XC, I  I , NOS 
BJLAL  NLQ1  Lo  i  t  E  (  lo)  ,  Hi  {  16)  ,  Pl(  Id  )  ,  IL(  16,27)  ,  oL  { Z  /)  ,  LI,  3,L  » 
1RR,C»Z,A,E, 

.  lTEMP,SUMi A 

_ LLU _ Lo _ I  =  I  t  N  E _ 

IF l ABSI X (R, I )) ,LE. 10 .0 )  GOTO  lo 

ftKi r e i 6,  i  n 

STOP 

16  C  0  N  r  I N  U  E 

I  f  FORMAT  I 'OSTATE  VARIABLE  BOUND  EXCEEDED*) 

_ 4=  X.1  R ,  1 ) _ 

A=u ( R  » I ) 

E  =  € I  R,2) 

6010(11,12)  1 1\« 

II  CO  13  1=1, NE 
G( I ) =0  •  C 

_ Cu  -  14  J  =  l,ivE _ 

14  J  X  I  I  ,  J  )  =  0 . 0 
DO  13  J=  1 ,  NK 

13  JK4I »J 1-0.0 
N  A=  2 

DO  15  i=l,NDS 

_ R E  AC  I  3 , 1  )  Fi  I)  ,N0(  1)  ,PI ill  ,PEII) _ 

wRiTE{6,i)Fm,Nom,pim,PE(n 
K  E AD  I  3  » 1  )  ITU  1, J)  *  J=  1 , 27 ) 

15  WRITE! 6* li (TL( I. J1 ,  J  =  1 , 2  7  ) 

1  FORMAT !8X,9F8*2) 

S  =  / •  4  L  2  4  2 

_ L=C.o9Cll6 _ 

RR=6236 1 CCO. C 
0=0,026543 
PL  4 1 J=0, 0 
DO  16  0=1*27 
CO  If  1=1, NDS 

_ 17  TL1I,J)=1TLU,J)  ->-460.0  )/  1800 .0 _ 

IP  {  J  •  E 0 •  1  )  GOTO  16 
DL I J ) =CL I J-l ) +0 
lo  OUNTIN0E 
12  J= 1,0+T/D 

T  E  M  P=  (  T  L  I  1 1 tJj^IDLlJ  +  ll-Ti+TLI  1I,J  +  1)*(T-DLIJ)  )  )/D 

_  P  =  F  I  I  I  1 ) -T*( PI  I  1  I )-PE(  I  1)  )/L _ 

SUMA=$*£XP  l  A-E/  TEMP  )  *P/  l  RR*TEMP*F(  I  I  )  ) 

G( i )=SUMA*( 1  ,C-Z)/ !i ,0+Nol I  I ) /F! I  I )+Z) 

IFl  I  JALK.LT.C)  GOTO  15 
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JX(  1»1)--$UMA*<2.  Q+NG  l  I  I  J/F(  I  I  )  )  /  (  1. 0+NU(  II) /H  (  I  I  )  +Z ) 
12 

Mi  1*U=G<1) 

JKl 1 »  2  ) =~G ( 1)/TEMP 

15  RETURN _ 

END 
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SUBROUTINE  SLOPE <R,T, I  JACK) 

C  SLOPE  SUBROUTINE  FOR  THE  PROPANE  PROBLEM  WITH  THE 

C  ORDER  OF  REACTION  UNSPECIFIED. 

REAL  X( iOflO ) » J  X  ( I0»10)»JK( 1 0 ♦ 1 0  > ,  G  { 10) ,  T  ,  C  ( II, 10) ,  X  0  { 3 
10,101 

I  f i  Ilj: S.R  i  N E-t  NK  ,  NA  ,  1  i  ,  NUS  ,  I  J  ACK _ 

COMMON  / AB 3/ C , N E , NK , NA/ AB5 / JX , JK , X ,G/AB8/X0, I  I  ,  NDS 
REAL  N0(  16) , F  (  16) , PI( 16)  ,PE(  16) ,TL ( 16,27) ,DL< 2  7 ) , D,S,L , 
1  RR  ,  Q  ,  Z  ,  A  ,  L  , 

I T  EMP , SUMA 
DO  16  I  =  1 , NE 

- IFiARS(X(RtI)).l  F  .  1  o  .  f)  )  GOTO  1A _ 

wR  I  T E  (  6 , 17) 

STOP 

16  CONTINUE 

17  FORMAT < •OSTATE  VARIABLE  BOUND  EXCEEDED* ) 

Z=X  (  R  ,  1  ) 

_ A=£_liUJa _ _ _ _ _ _ 

E=C(R, 2) 

G=C I R, 3 ) 

G'J  TO  {  1 1  ,  1 2  )  ,  NA 

11  DO  13  1=1, NE 
G ( I  1=0.0 

- DO  14  J=1,ME 

14  J  X (  I , J )  =  0 . 0 
DO  13  J= 1 , NK 

13  JK( I , J ) =0. a  _ 

N  A=  2 

DO  16  1=1, NOS 

_ R£A0(5,  1  1-  E(  I  )  ,NCi(  1  )  ,PI  (  I  )  ,PE(  I  )  _ 

WR I T  t ( 6 , 1  )  F {  I ) , NO (I),PI(I),PE(I) 

RE AO ( 5 , 1  )  ( T  L  < I  ,J)  ,  J  =  1 , 2  7  ) 

16  MR I T  E ( 6 , 1 ) ( T  L i I , J )  , J = 1 , 2  7 ) 

1  FORMAT ( 8 X,9F 8.2 ) 

S=7. 42242 

_ L=Q .6901 1 8 _ 

RR=62361 000 , 0 
0=0.026543 

DL (  1 ) =  0 . 0 

DO  16  J  =  1 , 2  7 
DO  17  1=1, NDS 

_ 17  TL(I,J)=(TL(  I,  J) +460. 0  1/1  800.0 _ 

IF(J.EQ.l)  GOTO  16 
DL ( J ) =DL( J-l  )+D 
16  CONTINUE 

12  J-l.O+T/D 

TEMP=<  TL  (  1 1  ,  J)  *.{  DL(  J  +  l  )-T  )+TL(  I  I  ,  J  +  l  )  *(  T-D  L  (  J  )  )  )  /D 
_ P=P I (II)- T*( PI  ( I  I ) -P  E ( I  I  )  )  / L _ 

SUMA-S*EXP< A-E/TEMP) /Fi I  I  ) 

G(  1  )  =  S  U  M  A  *  (  p  *  {  1.0-Z)  /(RR*TEMP*(  1 .0+NQ(  II)/F(II)+Z)))  **Q 


.  . 


(  ,  .  ,  )  -! 

,  j  ,  ,  i  )  ((ii  )  •  (  t  n  i. , . 

t(if 

,  ..... 

, ' ,  .  •  \  »  » 
t  .  }  .  ■  . .  )  . 

.  .  .  » 

■ 

i  .  .  .((.*..-■) 

(  ,  III 

'  ,  ! 

i  1  ,  )  = 

(  ,  )  -  ' 

(  ,  )  J  = 

(  »  ) 

,1 

.  =  .  )  ■  -1 
,  ,  ' 

,  [  -  i  .1 

(  i  ,t;)  ,  )  t  .I)  (!♦*) 

•  (  :  (i.  ) 

(  ,  ,  I  L  ,  )  <  r  C  ) 

(  .  f  -  .  ,  t  ,  ;  J  :  •  )  i  .  ;  l  M 

(  .  .  I  .  ' 

<  tv .  \  - 

(  i  •  C  --J 

. 

^  ,  I 

.  I  -  \  1  f 

\  (  .  •  (  L,  .  1  )  i  =  t  »  i  )  J  T  \  ' 

(  .  .  t  M  I 

M  -L  =  (. 


\  I  +  .  J  .. 

.  .  .  I  -  )  •  -  l  J 

-  (  1  ’  1  i  -  .  !  )  ! 

\  - 

H  \  l  I  /  »  .  '  )  \  1  ~  ,  )  «.  -  (  )  • 


020 


IF  (  I  JACR.LT.0)  GOTO  15 

J  x  C  1 , 1 )  =  -SUM A*Q* ( 2.0+NOI I  I )/F (  1 1 )  ) *( 1 . 0—  Z ) **{ Q- 1 .0)*  <P/ 
1  (ROTE HP  )  )  * 

1  *Q  /  (  1 . 0  +  N  i  J  (  II  )  /  F  (  I  I  )  +  z  )  *  *  (  1  .  o  +  Q  ) 

JK( 1 , 1 ) =G  C 1 ) 

- JN(  1,2)  =  -G  t  1 )  /  T  £  M  £ 

JK ( 1 »  3 ) -G ( 1 ) -ALOG !P*il .0-Z ) / < RR*TEMP*(  1 .0+NQ{  I  I)  /F ( I  I)  + 
1Z)  )  ) 

IS  RETURN 
END 
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APPENDIX  D 

Catalytic  Hydrogen  Chloride  Oxidation* 


The  data  obtained  for  the  oxidation  of  hydrogen 
chloride  was  taken  from  a  system  which  was  essentially 
a  batch  reactor.  Once  the  charge  entered  the  reactor, 
the  system  was  free  from  external  interference  except 
for  the  agitation  produced  by  a  circulating  pump  and  the 
removal  of  the  samples  for  analysis. 

An  adecruate  mathematical  model  of  this  system  is 
given  in  the  following  eciuations . 
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*  All  the  raw  data  for  this  example  was  supplied  by  Harding 


Bliss  of  Yale  University. 
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=  known  eauilibrium  constant 

=  partial  pressure  of  water 

=  partial  pressure  of  oxygen 

=  partial  pressure  of  chlorine 

=  partial  pressure  of  hydrogen  chloride 

=  total  number  of  gm-moles  initially  present  in  the 
reactor 

=  number  of  moles  of  oxygen  initially  in  the  reactor 
=  number  of  moles  of  HC1 
=  number  of  moles  of  chlorine 

=  number  of  moles  of  water  initially  in  the  reactor 
=  weight  of  catalyst  present  in  the  reactor. 

=  the  initial  pressure  of  the  reactor 
=  the  total  pressure  of  the  reactor  at  time  t 
=  number  of  moles  of  inerts 
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The  correction  factor  B(t)  was  correlated  from 
the  actual  pressure  data  using  polynomials  of  the  form- 

B  ( t )  =  bx  +  b2t  +  b3t2  +  b3t3  (D-6) 

The  maximum  amount  of  correction  was  about  fifteen  per 
cent  and  the  correction  factors  are  all  correct  to  within 
one  and  a  half  percent. 

All  the  data  used  to  identify  the  parameters  at 
355°C  are  shown  in  the  tables  below  starting  with  Table 
D-l .  The  data  used  to  identify  the  parameters  at  340°C 
are  shown  in  the  tables  begginine  with  Table  D-ll.  The 
data  used  to  attempt  the  identification  of  the  parameters 
at  325°C  are  shown  starting  with  Table  D-23. 


Data  for  Hydrogen  Chloride  Problem  at  355 
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TABLE  D-2 

Hydrogen  Chloride  Problem  Data  Set  One  at  355°C 


t  (min) 

P  (atm) 
u2 

P ( t ) (atm) 

0 

0.1501 

1.310 

1 

0.410 

1.272 

4 

0.118 

1.195 

7 

0.0958 

0.155 

10 

0.0853 

1.122 

13 

0.0772 

1.097 

TABLE  D-3 

Hydrogen 

Chloride  Problem  Data 

Set  Two  at  355 

t  (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.1496 

1.304 

1 

1.1404 

1.282 

4 

0.1150 

1.209 

7 

0.1002 

1.170 

10 

0.0890 

1.136 

13 

0.0795 

1.110 

21 

0.0633 

1.055 
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Hydrogen 

TABLE  D-4 

Chloride  Problem  Data 

Set 

Three  at  355 

t (min) 

P  (atm) 
u2 

P  (t)  (atm) 

0 

0.1495 

1.303 

1 

0.1380 

1.281 

4 

0.1155 

1.206 

7 

0.0998 

1.161 

Hydrogen 

TABLE  D-5 

Chloride  Problem  Data 

Set 

Four  at  355° 

t (min) 

P  (atm) 

P (t) (atm) 

0 

u2 

0.1520 

1.335 

1 

0.1445 

1.313 

4 

0.1147 

1.220 

7 

0.0988 

1.177 

10 

0.0890 

1.142 

13 

0.0808 

1.119 
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TABLE  D-6 


Hydrogen  Chloride 

Problem  Data 

Set  Five  at  335°' 

t (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.1560 

1.322 

1 

0.4196 

1.298 

4 

0.1252 

1.221 

7 

0.1115 

1.183 

10 

0.1006 

1.152 

13 

0.0907 

1.124 

19 

0.0785 

1.080 

TABLE  D-7 

Hydrogen  Chloride 

Problem  Data 

Set  Six  at  35 

t (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.2221 

1.330 

1 

0.2095 

1.310 

4 

0.1831 

1.250 

7 

0.1720 

1.226 

10 

0.1660 

1.210 

13 

0.1622 

1.201 
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TABLE  D-8 


Hydrogen 

Chloride 

Problem  Data 

Set 

Seven  At  355 

t  (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.2160 

1.311 

1 

0.2050 

1.305 

4 

0.1809 

1.245 

TABLE  D-9 

Hydrogen 

Chloride 

Problem  Data 

Set 

Eight  at  355 

t  (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.0745 

1.330 

10 

0.0491 

1.234 

TABLE  D-10 

Hydrogen 

Chloride 

Problem  Data 

Set 

Nine  at  355° 

t  (min) 

P  (atm) 
u2 

P ( t) (atm) 

0 

0.1405 

1.317 

2 

0.1225 

1.270 

5 

0.1060 

1.211 

8 

0.0949 

1.170 

11 

0.0882 

1.140 

14 

0.0788 

1.110 

21 

0.0654 

1.054 

i  . 

Data  for  Hydrogen  Chloride  Problem  at  340 
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TABLE  D-ll  (continued) 
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TABLE  D-12 

Hydrogen 

Chloride 

Problem  Data 

Set 

One  at  340°c 

t  (min) 

Pn  (atm) 
u2 

P ( t) (atm) 

0 

0.1500 

1.307 

1 

0.1416 

1.284 

4 

0.1260 

1.217 

7 

0.1123 

1.182 

13 

0.0967 

1.133 

TABLE  D-13 

Hydrogen 

Chloride 

Problem  Data 

Set 

Two  at  340°C 

t  (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.1512 

1.319 

1 

0.1418 

1.297 

4 

0.1228 

1.234 

7 

0.1147 

1.195 

TABLE  D-14 

Hydrogen 

Chloride 

Problem  Data 

Set 

Three  at  340 

t (min) 

P  (atm) 

U2 

P ( t) (atm) 

0 

0.1530 

1.334 

1 

0.1426 

1.321 

4 

0.1269 

1.266 

7 

0.1126 

1.266 

13 

0.0988 

1.102 

• 
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TABLE  D-15 


Hydrogen  Chloride 

Problem  Data 

Set 

Four  at  340 

t  (min) 

P  (atm) 
u2 

P  (t)  (atm) 

0 

0.1535 

1.340 

1 

0.1475 

1.322 

4 

0.1241 

1.253 

7 

0.1140 

1.216 

TABLE  D-16 

Hydrogen  Chloride 

Problem  Data 

Set 

Five  at  340 

t  (min) 

P  (atm) 

°2 

P ( t) (atm) 

0 

0.1548 

1.312 

1 

0.1482 

1.291 

4 

0.1329 

1.241 

13 

0.1078 

1.172 

22 

0.0925 

1.122 

TABLE  D-17 

Hydrogen  Chloride 

Problem  Data 

Set 

Six  at  340° 

t (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.1560 

1.322 

1 

0.1450 

1.310 

4 

0.1290 

1.245 

13 

0.1062 

1.165 

• 

-  D-13  - 

TABLE  D-18 


Hydrogen 

Chloride 

Problem  Data 

Set 

Seven  at  340 

t  (min) 

P  (atm) 

2 

P ( t) (atm) 

0 

0.1770 

1.301 

1 

0.1705 

1.281 

4 

0.1506 

1.228 

7 

0.1365 

1.192 

16 

0.1123 

1.114 

TABLE  D-19 

Hydrogen 

Chloride 

Problem  Data 

Set 

Eight  at  340 

t  (min) 

P  (atm) 
u2 

P ( t) (atm) 

0 

0.2210 

1.320 

1 

0.2155 

1.310 

7 

0.1825 

0.232 

13 

0.1691 

1.200 

TABLE  D-20 

Hydrogen 

Chloride 

Problem  Data 

Set 

Nine  at  34  0°C 

t  (min) 

P  (atm) 

U2 

P ( t) (atm) 

0 

0.2205 

1.319 

1 

0.2085 

1.295 

7 

0.1762 

1.215 

13 

0.1662 

1.193 
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TABLE  D-21 

Hydrogen  Chloride  Problem  Data  Set  Ten  at  340°C 

t(min)  P  (atm)  P(t)(atm) 

U2 

0  0.0843  1.316 

10  0.0694  1.225 

TABLE  D-2 2 

Hydrogen  Chloride  Problem  Data  Set  Eleven  at  340°C 

t(min)  Pn  (atm)  P(t) (atm) 

u2 

0  0.1535  1.339 

1  0.1475  1.322 
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Hydrogen 

TABLE  D-24 

Chloride  Problem  Data 

Set  One  at  325 

t  (min) 

P  (atm) 
u2 

P ( t) (atm) 

0 

0.1560 

1.333 

1 

0.1497 

1.311 

4 

0.1392 

1.255 

7 

0.1335 

1.222 

10 

0.1271 

1.191 

13 

0.1238 

1.167 

Hydrogen 

TABLE  D-25 

Chloride  Problem  Data 

Set  Two  at  325 

t (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

Q . 1566 

1.327 

1 

0.1503 

1.295 

4 

0.1380 

1.252 

7 

0.1338 

1.236 

10 

0.1277 

1.220 

13 

0.1238 

1.210 

18 

0.1181 

1.191 

D-17 


TABLE  D-26 

Hydrogen  Chloride  Problem  Data  Set  Three  at  325°C 


t (min) 

P  (atm) 
u2 

P ( t) (atm) 

0 

0.2220 

1.344 

2 

0.2020 

1.317 

5 

0.1840 

1.263 

10 

0.1746 

1.251 

15 

0.1674 

1.237 

20 

0.1631 

1.227 

TABLE  D-27 

Hydrogen 

Chloride  Problem  Data 

Set  Four  at  325° 

t (min) 

P  (atm) 
u2 

P ( t) (atm) 

0 

0.2180 

1.320 

5 

0.1850 

1.250 

10 

0.1719 

1.222 

15 

0.1654 

1.208 

20 

0.1604 

1.195 

D-18 


TABLE  D-28 

Hydrogen  Chloride  Problem  Data  Set  Five  at  325°C 


t (min) 

P  (atm) 

2 

P (t) (atm) 

0 

0.0748 

1.333 

1 

0.0741 

1 .327 

4 

0.0706 

1.306 

7 

0.0710 

1.298 

10 

0.0701 

1.293 

13 

0.0694 

1.280 

20 

0.0648 

1.246 

Hydrogen 

TABLE  D-29 

Chloride  Problem  Data 

Set  Six  at  325° 

t  (min) 

P  (atm) 
u2 

P (t) (atm) 

0 

0.1417 

1.328 

2 

0.1362 

1.298 

5 

0.1268 

1.255 

8 

0.1193 

1.227 

11 

0.1160 

1.197 

14 

0.1130 

1.170 

23 

0.1005 

1.106 

D-19 


TABLE  D-30 

Hydrogen  Chloride  Problem  Data  Set  Seven  at  325°C 


t  (min) 

Pn  (atm) 
u2 

P (t) (atm) 

0 

0.1415 

1.320 

1 

0.1378 

1.297 

4 

0.1228 

1.240 

7 

0.1165 

1.210 

10 

0.1101 

1.189 

13 

0.1061 

1.167 

19 

0.0989 

1.132 

TABLE  D-31 

Hydrogen  Chloride 

Problem  Data 

Set  Eight  at  325 

t  (min) 

Pn  (atm) 
u2 

P (t) (atm) 

0 

0,1420 

1.326 

1 

0,1380 

1.304 

4 

0.1243 

1.255 

7 

0.1177 

1.228 

10 

0,1112 

1.204 

13 

0,1056 

1.186 

19 

0.0981 

1.150 

D-2  0 


TABLE  D-32 


Hydrogen  Chloride  Problem 


t(min)  P  (atm) 

U2 


0 

0.1410 

1 

0.1390 

4 

0.1235 

7 

0.1122 

10 

0.1120 

13 

0.1064 

Data  Set  Nine  at  325°C 

P (t) (atm) 

1.319 

1.304 

1.256 

1.235 

1.221 

1.206 


An  attempt  was  made  to  determine  the  unknown 
parameters  in  equation  (D-l)  at  each  of  the  three  tem¬ 
peratures.  The  results  of  each  of  these  attempts  is 
presented  below.  At  355°C  the  solution  coverged  direct¬ 
ly  with  quasilinearization  as  shown  in  Table  D-33. 
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TABLE  D-3 3 

Solution  of  Hydrogen  Chloride  Problem  at  355°C 


Iteration 

Number 

a^xlO^ 

a2 

a3 

a4 

Sum  of 
Errors 
Squared 

0 

2.80 

3.10 

25.0 

230 

0.00315 

1 

5.44 

6.91 

36.1 

401 

2 

9.00 

12.2 

50.0 

661 

3 

9.47 

12.6 

47.3 

704 

4 

9.55 

12.8 

47.0 

732 

5 

9.76 

13.1 

47.8 

755 

6 

9.78 

13.1 

47.8 

759 

7 

9.80 

13.2 

47.9 

762 

8 

9.80 

13.2 

47.9 

762 

0.000747 

At 

34 0°C  a 

solution 

could  also 

be  obtained  using 

quasilinearization  directly  as  shown  in  Table  D-34. 
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TABLE  D-34 

Solution  to  Hydrogen  Chloride  Problem  at  340°C 


Iteration 

Number 

a^xlO^ 

a2 

a3 

a4 

Sum  of 
Errors 
Squared 

0 

0.920 

0.220 

21.0 

142 

0.00708 

1 

2.79 

0.651 

90.3 

-65.1 

2 

5.28 

0.732 

128.8 

-50.3 

3 

5.57 

0.793 

120.6 

-34.8 

4 

5.58 

0.793 

121.6 

-34.9 

5 

5.57 

0.790 

121.3 

-34.8 

6 

5.57 

0.791 

121.3 

-34.8 

0.00144 

At  325°C  no  solution  was  found.  An  attempt  to 

find  the  solution  directly  using  quasilinearization  failed 

because  the  maximum  allowable  relative  change  in  the 

parameter  C  was  exceeded.  The  initial  values  of  S 

max  max 

and  C  were  0.05  and  50  respectively.  The  results  are 
max  ^  J 

summarized  in  Table  D-35.  The  SLOPE  subroutine  used  for 
this  problem  is  listed  at  the  end  of  this  section. 
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APPENDIX  E 


Chemcell  Problem 


Batch  reactor  data  has  been  accumulated  for  a 
complex  kinetic  system.  Six  experiments  were  performed 


at  100°C,  and  five  at  80°C.  Three  principal  components 


x^ ,  x^/  and  x^  can  be,  and  were  directly  measured.  A 
fourth  is  given  by  a  linear  relation  between  x^ ,  x^  and 
one  or  more  intermediates  which  exist  but  cannot  be  measured. 

Each  of  the  experiments  was  performed  in  a  dilute 
acidic  solution.  The  concentration  of  water  (W)  and  of 
hydrogen  ion  (H+)  were  constant  throughout  any  one  experi¬ 
ment.  The  initial  conditions  for  each  of  the  components 
except  x^  were  all  zero.  The  initial  condition  for  x^  was 
the  constant  c  . 


o 


The  first  model  that  was  proposed  could  be  repre¬ 
sented  by  the  system  of  eauations  given  below. 


dx 


1 


dt 


(E-l) 


dx 


2 


a2^x4  a3X5 


(E-2) 


dt 


E-2 


dx. 

dt 


a4X5 


(E-3) 


+  +  + 

(a7Hx0  -  a0W)a1Hx1  -  a-,Hx0  (c  -  x.,  -  x^)  (a^  +  a0W) 


„  _  '7  2  -2"'~1“~1  “7“"2  'w0 

5  ~  - T -  - 


(a^W  -  a7Hx2)a6Xg  -  (a^x2  +  a^  +  +  a7Hx2) (a,.  +  a2W) 


x.  + 
4 


+ 

a^Hx^ 
a5  +  a2W 


a6X2 


ac  +  a„W 
5  2 


x. 


(E-4 ) 
(E-5) 


x6  =  Co  “  (X1  +  x3  +  X4  +  X5} 


(E-6) 


All  attempts  to  fit  the  above  model  failed  to  con¬ 
verge  . 

A  second  and  partially  successful  attempt  was  made 
to  fit  the  data  to  the  model  given  below. 


dx 


dt 


-  =  -  a1H+Wx1 


(E-7 ) 


dx 


dt 


-  =  a2x4  ‘  a5H+x2x6  +  alH+Xl 


(E-8 ) 


dx. 


dt 


a3X4 


(E-9) 


x 


a.x.H  W  +  acx0(c  -  x  -  x_.) 
1  1_ 5  2  o  1  o 

a~  ■+•  a t  a [- H  x0 
2  J  5  z 


(E-10 ) 


E-3 


A  different  solution  was  obtained  for  this  prob¬ 
lem  for  each  initial  guess  employed.  With  the  80°C 
data,  different  solutions  were  found  for  each  initial 
guess.  The  first  solution  was  found  starting  with  the 
initial  guess: 


al 

“0.10  ~ 

a2 

1.00 

a3 

0.10 

_a5_ 

1.00 

and  with  the  initial  maximum  step  size  Smax  being  0.01. 
With  this  initial  guess  the  results  in  Table  E-l  were 


obtained . 


. 


Phas 

Step 

0 

2 

3 

4 

5 

6 

7 

8 

9 

Phas 

Iter 

No 

0 

1 

2 

3 

4 


E-4 


TABLE  E-l 


First  Solution  for  80°C  Data 


Sum  cf 


al 

a2 

a 

3 

a5 

Errors 

Scruared 

0.10 

1.00 

0.10 

1.00 

0.3475 

0.0848 

0.976 

0.131 

0.508 

0 . 3063 

0.0715 

0.956 

0.157 

0 . 245 

0 . 2660 

0.0701 

0.940 

0.179 

0.0919 

0.2266 

0.0504 

0.929 

0.196 

0.00127 

0.1889 

0 . 1421 

0.923 

0.209 

-0.0474 

0.1528 

0.0351 

0.920 

0.217 

-0.0626 

0.1187 

0.0291 

0.923 

0.220 

-0.0453 

0.08703 

0.0161 

0.963 

0.191 

0.335 

0.01584 

0.0161 

0.903 

0.191 

0.335 

0 .00963 

1.100 

0.154 

0.822 

0.0106 

1.019 

0.136 

1.249 

0.01064 

1.017 

0.1379 

1.302 

0.01064 

1.017 

0.1379 

1.301 

0.00007009 


• 

T3 

E-5 


With  the  second  initial  guess 


0.040 

1.90 

0.58 

1.58 


(E-12 ) 


and  the  same  value  for  the  results  were  as  presented 

max  r 

in  Table  E-2. 


TABLE  E-2 

Second  Solution  for  80°C  Data 


Phase  I 

Step 

No. 

al 

a2 

a3 

a4 

Sum  of 
Errors 
Sauared 

0 

0.040 

1.90 

0.580 

1.580 

0.1351 

1 

0.0329 

1.896 

0.600 

1.120 

0.1011 

2 

0.0175 

2.007 

0.516 

0.786 

0.02112 

Phase  II 
Iteration 
No. 

0 

0.0175 

2.007 

0.516 

0.786 

1 

0.00908 

2.149 

0.381 

0.745 

2 

0.0107 

2.244 

0.284 

1.189 

3 

0.0107 

2.227 

0.302 

1.298 

4 

0.0106 

2.227 

0.302 

1.293 

0.00006993 

- 

. 

E-6 


With  the  third  initial  guess 


0.1 


a 


0.1 


0.1 


(E-13) 


0.1 


and  with  the  same  initial  value  for  S  the  results  in 

max 


Table  E-3  were  obtained. 


TABLE  E-3 

Third  Solution  for  80°C  Data 


Phase  I 
Step  No. 

al 

a2 

a3 

a4 

Sum  of 
Errors 
Sauared 

0 

0.1 

0.1 

0.1 

0.1 

0.2553 

1 

0.0766 

0.102 

0.0996 

-0.0282 

0.2154 

2 

0.0585 

0.106 

0.0978 

-0.1098 

0.1740 

3 

0.0448 

0.1106 

0.0944 

-0.1547 

0.1332 

4 

0.0345 

0.1173 

0.0891 

-0.1678 

0.09506 

5 

0.0154 

0.01547 

0.0535 

0.1374 

0.01070 

Phase  II 
Iteration 
No. 

0 

0.0154 

0.1547 

0.0535 

0.137 

1 

0.00979 

0.1780 

0.0305 

0.561 

2 

0.0106 

0.1840 

0.0245 

1.120 

3 

0.0106 

0.1835 

0.0245 

1.377 

4 

0.0106 

0.1835 

0.0250 

1.374 

0.00007195 

E-7 


A  similar  phenomenon  was  encountered  when  attempts 
were  made  to  determine  the  parameters  characteristic  of 
the  100°C  data.  Three  attempts  were  made,  all  with 
different  results.  For  each  of  the  attempts  the  initial 
value  of  the  maximum  step  size  was  retained  at  0.01. 

For  the  first,  second,  and  third  initial  guesses, 
the  results  are  given  in  Tables  E-4,  E-5,  and  E-6  respec¬ 
tively  . 

TABLE  E-4 

First  Solution  for  100°C  Data 

Sum  of 


Step  No. 

al 

a2 

a3 

a4 

Errors 

Squared 

0 

0.1 

1.0 

0.1 

1.0 

0.2470 

1 

0.0124 

1.03 

0.0760 

17.7 

0.2048 

2 

0.0152 

1.04 

0.0619 

33.6 

0.1664 

3 

0.0185 

1.05 

0.0582 

41.5 

1.1318 

4 

0.0223 

1.05 

0.0644 

38.0 

0.1009 

5 

0.0423 

1.00 

0.121 

13.9 

0 .03033 

6 

0.0692 

0.975 

0.151 

7.61 

0.001763 

7 

0.0815 

0.979 

0.150 

7.14 

0.0002298 

’ 

.. 

E-8 


TABLE  E-5 

Second  Solution  for  100°C  Data 


Step  No. 

al 

a2 

a3 

a4 

Sum  of 
Errors 
Scruared 

0 

0.070 

2.00 

0.400 

5.00 

0.002076 

1 

0.0787 

2.06 

0.339 

6.52 

0.0003089 

2 

0.0815 

2.10 

0.321 

7.11 

0.0002301 

TABLE  E-6 

Third 

Solution  for  10Q°G  Data 

Step  No. 

ai 

a2 

a3 

a4 

Sum  of 
Errors 
Swuared 

0 

0.1 

0.1 

0.1 

0,1 

0,02516 

1 

0.0922 

0.126 

0.0738 

1,30 

0,009845 

2 

0.0854 

0.154 

0.0461 

3.12 

0,001711 

3 

0.0815 

0.174 

0.0270 

7.47 

0.0002271 

The  data  used  t©  produced  the  results  in  the 
previous  si*  tables  are  listed  in  Tables  E^7  through  to 
E-17 ,  The  SLOPE  subroutines  used  are  listed  immediately 
following  the  data  * 


E-9 


TABLE  E-7 


Data  Set  One 

u 

o 

o 

00 

(H+)  =  1.27  x 

10  ^  g-mole/1 

c  =  7.91  x 

o 

10  ^  g-mole/1 

t  (min) 

X1 

X2  9 

X- 

X 

(g-mole/1) xl0z 

j 

0.0 

7.91 

0.00 

0.00 

5.3 

7.61 

0.52 

0.07 

9.8 

7.39 

0.97 

0.07 

19.9 

6.87 

1.93 

0.14 

29.8 

6.40 

2.87 

0.14 

43.7 

5.79 

4.04 

0.20 

59.6 

5.08 

5.31 

0.34 

80.1 

4.41 

6.59 

0.41 

100.5 

3.78 

7.78 

0.48 

150.3 

2.68 

9.85 

0.62 

204.9 

1.76 

11.46 

0.83 

14440.0* 

0.00 

13.07 

2.75 

28880.0* 

0.00 

11.81 

4.01 

*These  points  were  not  used  to  find  the  parameters 


1  . 

E-10 


TABLE  E-8 
Data  Set  Two  80°C 


(H+)  = 

2.30 

x  10  2  g-mole/1 

c  = 

o 

7.95 

x  10  2  g-mole/1 

t  (min) 

X1 

x2 

(g-mole/l)xl02 

X3 

0.0 

7.95 

0.00 

0.00 

2.57 

7.69 

0.44 

0.07 

5.32 

7.40 

1.03 

0.07 

10.13 

6.92 

1.98 

0.07 

15.20 

6.49 

2.77 

0.14 

20.05 

6.06 

3.60 

0.19 

29.70 

5.26 

5.12 

0 . 27 

39.82 

4.55 

6.38 

0.41 

59 . 50 

3.47 

8.41 

0.53 

80.00 

2.64 

9.94 

0.67 

101.47 

1.96 

11.17 

0.80 

150.65 

0.98 

12.91 

1.01 

200.48 

0.50 

13.67 

1.23 

14440.0* 

0.00 

11.98 

3.92 

28880.0* 

0.00 

10.66 

5.23 

* 


These  points  were  not  used  to  find  the  parameters 


. 

, 

. 

E-ll 


TABLE  E-9 

Data  Set  Three  at  80°C 


<H+) 

=  3.68 

x  10  2 

g-mole/1 

c 

o 

=  8.22 

x  10~2 

g-mole/1 

t  (min) 

X1 

x2  o 

x. 

JL 

<g- 

-mole/1 ) xlO 

0.0 

8.22 

0.0 

0.00 

1.75 

7.93 

0.51 

0.07 

4.72 

7.48 

1.41 

0.07 

10.52 

6.55 

3.13 

0 . 20 

14.87 

6.02 

4.19 

0.20 

19.78 

5.37 

5.35 

0.34 

30.28 

4 . 32 

7.34 

0.47 

39.82 

3.4  8 

8.87 

0.60 

61.00 

2.20 

11.23 

0.80 

82.12 

1.35 

12.79 

0.96 

99.98 

0.92 

13.50 

1.09 

149.95 

0.32 

14.38 

1.43 

299.40 

0.11 

14.57 

1.64 

14440.00* 

0.00 

11.28 

5.15 

28880.00* 

0.00 

10.12 

6 . 32 

*  These  points  were  not  used  to  find  the  parameters 


E-12 


TABLE  E-12 

Data  Set  ^our  at  80°C 


(H+)  =  3.45 

x  10  2  g-mole/1 

CQ  =8.25 

x  10  2  g-mole/1 

t  (min) 

X1 

X2 

(g-mole/1) xlO2 

X3 

0.0 

8.25 

0.00 

0.00 

2.42 

7.86 

0.73 

0.07 

4.85 

7.49 

1.38 

0.14 

9.97 

6.71 

2.88 

0.20 

15.03 

6.00 

4.18 

0.33 

19.80 

5.44 

5.28 

0.34 

29.63 

4.44 

7.10 

0 . 53 

41.43 

3.42 

8  99 

0.68 

59.55 

2.38 

10.95 

0.81 

83.95 

1.41 

12.66 

1.04 

101.30 

0.96 

13.45 

1.14 

150.55 

0.39 

14.23 

1.50 

202.75 

0.14 

14.50 

1.72 

14440.00* 

0.0 

11.31 

5.19 

28880.00* 

0.0 

9.51 

6 . 99 

* 


These  points  were  not  used  to  determine  the  parameter 


• 

' 

• 

E-13 


TABLE  E-li 


Data 

Set  Five 

at  80°C 

(H+) 

=  6.98  x 

10  2  g-mole/1 

°o 

=  8.35  x 

10  2  g-mole/1 

t  (min) 

X1 

g-mole^l) xlQ2 

x3 

0.00 

8.35 

0.00 

0.00 

7.40 

7.60 

1.35 

0.14 

9.80 

5.64 

5.08 

0.34 

15.2 

4.45 

7.32 

0.49 

19.8 

3.71 

8.69 

0.60 

24.9 

3.01 

4.99 

0.70 

29.9 

2.41 

11.08 

0.80 

40.0 

1.63 

12.49 

0.95 

50.1 

1.05 

13.41 

1.18 

60.2 

0.67 

14.06 

1.79 

87.5 

0.22 

14.71 

1.55 

104.7 

0.10 

14.63 

1.86 

156.7 

0.04 

14.36 

2.27 

211.8 

0.00 

13.95 

2.76 

14440.0 

0.00 

10.26 

6.44 

28880.0 

0.00 

9.73 

6.97 

. 

E-14 


TABLE  E-12 


Data  Set  One  at  100°C 


(H+) 

=  6. 

14  x 

10  3  g-mole/1 

c 

o 

=  8. 

28  x 

10  ^  g-mole/1 

t  (min) 

X1 

*2 

x 

(g-mole/1 ) xlO 

j 

0 

8.28 

0.00 

0.00 

5.9 

7.08 

2.26 

0.14 

10.6 

6.17 

3.95 

0.27 

40.1 

2.59 

10.49 

0.89 

59.9 

1.44 

12.66 

1.03 

100.2 

0.44 

14.30 

1.37 

149.9 

0.11 

14.68 

1.66 

206.0 

0.03 

14.44 

2.06 

TABLE 

E-13 

Data 

.  Set 

Two 

at  100°C 

<H+) 

=  9 

.86 

x  10  3  g-mole/1 

c 

o 

=  7 

.81 

-2 

x  10  g-mole/1 

t  (min) 

X1 

x2  2 

X3 

(g-mole/1 ) xlO 

0.00 

7.84 

0.00 

0.00 

5.1 

6 . 28 

2.92 

0.20 

10.2 

4.88 

5.50 

0.42 

20.7 

2.91 

9.16 

0.70 

39.6 

1.16 

12.33 

1.04 

60.1 

0.40 

13.57 

1.30 

100.3 

0.07 

13.88 

1.66 

150.1 

0.00 

13.62 

2.06 

200.0 

0.00 

3  3.26 

2.43 

E-15 


t  (min) 

0.0 

5.1 

20.2 

39.8 
60.0 

99.8 

152.8 
201.7 


t  (min) 

0.0 

5.9 

20.5 

40.6 

59.7 

100.0 

150.8 

203.9 


TABLE  E-14 


Data  Set  Three  at  100°C 


(H+)  =  1.25 

x  10  2  g-mole/1 
-2 

c  =7.49 

o 

x  10  g-mole/1 

X1 

x2 

(g-mole/1) xlO2 

X3 

7.49 

0.00 

0.00 

5.80 

3.16 

0.21 

2.36 

4.57 

0.69 

0.75 

12.38 

1.09 

0.21 

13.28 

1.30 

0.04 

13.21 

1.73 

0.00 

12.80 

2.18 

0.00 

12.42 

2.56 

TABLE  E-15 

Data  Set  Four  at  100°C 

(H+)  =  2.41 

_3 

x  10  g-mole/1 

o 

0 

ii 

• 

4^ 

-2 

x  10  g-mole/1 

X1 

(g-mole^l ) xlO2 

X3 

7.74 

0.00 

0.00 

7.29 

0.82 

0.07 

6.23 

2.82 

0 . 20 

5.05 

5.05 

0.34 

4.06 

6.88 

0.48 

2.72 

9.35 

0.69 

1.60 

11.40 

0.88 

0.85 

12.69 

1.09 

■ 

« 

' 

e  .  i 

E-16 


TABLE  E-16 


Data 

Set  Four 

at  100°C 

(H+) 

=  1.03  x 

10  2  g-mole/1 

c 

o 

=  1.563  ; 

x  10  2  g-mole/1 

t  (min) 

X1 

x2 

(g-mole/1) xlO2 

X3 

0.0 

1.563 

0.000 

0.000 

3.1 

1.374 

0.353 

0.026 

6.0 

1.187 

0.699 

0.053 

15.5 

0.748 

1.514 

0.117 

24.7 

0.476 

2.014 

0.161 

41.3 

0.206 

0.474 

0.241 

62.3 

0.077 

2.642 

0.330 

88.0 

0.020 

2.670 

0.415 

108.1 

0.013 

2.628 

0.468 

157.3 

0.003 

2.528 

0.595 

210.0 

0.000 

2.431 

0.695 

E-17 


TABLE  E-17 

Data  Set  Six  at  100°C 


(H+)  = 

c  = 

o 

1.20 

7.59 

x  10  2  g-mole/1 
x  10  2  g-mole/1 

t (min) 

X1 

x2 

(g-mole/1) xlO2 

X3 

0.0 

7.59 

0.00 

0.00 

2.4 

6.84 

1.44 

0.07 

5.3 

5.92 

3.12 

0 . 22 

9.8 

4.59 

5.59 

0.41 

15.1 

3.47 

7.70 

0.54 

20.0 

2.64 

9.22 

0.68 

30.6 

1.47 

11.41 

0.83 

39.6 

0.90 

12.35 

1.04 

49.5 

0.52 

12.98 

1.17 

60.9 

0.30 

13.17 

1.40 

80.5 

0.11 

13.42 

1.54 

101.6 

0.00 

13.35 

1.83 

160.7 

0.00 

13.02 

2.16 
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5  0  c  k  ^  u  1  1 1\  t  $LCPb(K,T  ,  I  J  A  C  K  i 

,_X - S.L fl£g  A  J o  k  l U T  i  !\  L  £ Q *  T  n C  S£V£  4  CONSTANT  CHEMCELL 

^  PROBLEM  A  i  80  DEGREES* 

t<b  A  L  *4  X(10*10)fJX{lCtl0JtJK(10fl0)tG<10),T*C(ll»10)#XG 

_ I  (30, 10)  _ _ _ 

Il\l  LvjcK  Kjjsic  ,  N  K  »  N  Li  ,  1  1 ,  IM  0  S  ,  1 J  A  C  R 

CCrtMUN  /  A63  /  0  ,  iME,  (MK  ,  (MU/AB5/ JX  ,  JK  »  A  ,G/AbB/  XC,  1  i  ,  NUS 
R  £  a  L  '<r  -+  X4 1 X5 1 X6 1 K 1 «  K2  »K3«  K4 » K5  «  K6  .K7.HIQ5)  ,  W  ( Q 

154  »CQ{ 05 ) 

i  tt  F  , Gb  , Gu2  »  U  X4 X 1 , D x4  a2  »  C  X4 X  3 , D X4R 1 , DX4r2 ,D  X4K3  »  0X4R4 , DX 

_ l4R3,JX4Ko, _ _ _ ___ _ _ 

ILSX4K7  ,0X3X1  ,  Dxi>X2,CX5X3tDX5Kl,LX3K2,DX5K3,DX3K4,DX5K5,D 

1 X3RO , UX5K  f 
INTEGER  x  ,  j 
00  16  l ,  N  E 

IF  i  AdA ( X { R ,  L  1  )  ,Lc  .  10 .0  )  GO f U  16 

_ XKITlU,  1  /} _ _ 

6 1  C  P 

16  CUNT i N  L  b 

— _ U . EORMAIi AOSTA  IE  VARIABLE  BfiUiiC  EXCEEDED »  ) _ _ 

X  i  —  A  (  K  ,  1  ) 

X  2  =  X  l K  t  2  I 

_ X3=X(R,3J _ 

K  1  =  0  i  R  ,  1 ) 

K  2  =  0  I  R  ,  2  1 

& 3  =0  (R  *3) _ 

K4  =  0 IR  ,4  ) 

K  3  =  0  (  R  ,  5  ) 

_ R  o  =  0 ( R , 6  )  _ 

K  7  =  0  t  R  » 7  ) 

GOTH  11 ,12)  j  (\i, 

. il.Liu  13  1=1,  HE  _  , _ 

G  (  1 ) =0  .0 
UU  1  4  J  =  1  ,  I'M  E 

_ 1 4  J  X  (  1  i  J  )  =  u  «  J _ 

Lu  13  J  =  1  ,  N  K 
13  jM  I  ,  J  )  =0.0 
f\U=  A 

U  m  T  a  nil  )  , ill  tCu(  1 )  /0  •  0 1 27  *  55  •  6 , 0  •  0  79 1  / 

LATA  h  ( 2  J  t  >i  12  J  ,0L(  2  J  /  J  .  0  2  3  0  » 53.6,0.0793/ 

_ UAl A  ri  ( 3  >  ,  W  i  3  1  ,0c(3)/Q. 0  36 8,53 .6, 0,082  2/ _ 

U  A  T  a>  M  (  4  J  i  w  (  4  1  ,  C  C  l  4  J  /  0 .0343, 55 .6 ,0.0625/ 

UATm  H  (  5  1  ,  M ( 5  J  »C0{ 3J/0. 0696, 55.6,0.0635/ 

1 4.  PF=  (  K  /  *  H  l  III  *X2~K2*N1  113  )#K1*H(  II  )»X1~K  7*H  ill  >*X2*CCOC  1 
i  i  )  -  i  a  1  +  X  3  ;  ) 

1* ( K3+K2*X  i  1  I  )  1 

c  b-=  (  a  2  *  iM  (  11  )  ~  N  7  ^  H  {  I  1  1  *X2  )  *K6*x2-  iR6*X2+R3+R4+K  7*H<  I  i  1*X 
1 2  i  *  l  K5  +  K2# 

lrt(  1  1  I  i 
X3=FF/ bb 
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X4  =  (  k  1  *H  (  I  i  j^Xl+Kb^XZ^Xii/  (Kd-HvZ^kIII)  J 
S  {  l.  )  -  “  & l!  *  fcjj  i  i  i  )  x  i  £  K  5  *  X  ± 

G(  2  )^&2 III) *X4+K3*X5-K. 6 *X2*X5-K7*HC 1 1 1 *X2*(CG< II J- ( X 1 
1  +X3+X4+Xb  J  ) 

G  (  3  )  —  K  4  ^  X  3 _ _ _ _ _ _ 

IF(  IJACk.LT  «0)  bufU  15 

b/AjXi=  i  U7*rH  1 1  I  *  X  2  -  i\  2  *  «i  (  1  I  i  J  *  K 1  #  H  (  II  I  +  K  7  *  hi  (  II  )  *  X  2  *  (  K  5  +• 

i  X  2  *  (  1 1  i  1  )  /  _ _ 

IGG 

J  a  3  X  —  (  l  K  /  i\  I  X  1  ^  H  {  I  i  }  2  —  K  7  {  i  1  )*(Lu(  i  1  )-(  X  I  X  J  )  )  *  (  K  3 

L  +t\Z  )  */i  {  L  l  i  I _ _____ _ 

i  *Go  -  (KG  ’•‘  l  K2*  v\  l  i  1  )  -  kY  =rH  {  I  1  )  #X2  J  -K  7^H  (  I  I  j*K63i:X2--(K3  +  K2*W( 
liiJ  I  # ( K  5  +  K  7 
2*HI  i  i  )  U*FFl/GG**2 
DX5X3^K7*H{ 1 I J*X2*(K 5 +K2*W< I U )/GG 
JX3  K  i  =  I K  7  {  II  J - K2  * W (  I  1)  J  *  H  I  I  I  ) *  X 1 / uG 
_ QX-3K2-  ( ( -*( I  1  ) *k1*h ( 1  I  )»X  1-n (  1  1  )  »K7*h( 1  1  ) *X2»{ CGI  1  I ) -( X 

1 1 +X3  )  )  ) *oo 

1  -  4  W  I  1 1  J  *  i\  b  *  X  2  -  w  (  II  )*(K6*X2i-K3-*-K4  +  K7*h(  II  )  *  X  2  J  I^FFJ/GG*# 

JLZ. _ _ _ _ _ _ _ _ _ _ _ 

0  X  3  K  J  -  r  F { K  5  +  K  2  v  h  (  1 1  I  |  /GG  ♦ 2 
QX5K4=L:X3K3 

_ li-A-5.K3  =  i  -K/*r1  (Il)*X2*lLC(  1  U - <X 1  +  X3  )  )  *G^  +  I  Ku*X2  +  K  7 

i  ;<4  Fi  {  111^X21* 
iFFi/GG^rSZ 

u^jmxjh  KjLtK2*M  U  II  I-(K2*W1 1  I  )-K7*ii(  I1)*a2I  *X2I  *Fi  /G 
I G :'c  ii 

0X5K  7=1  (  X2*X1*K1*H  l  I  I)  **2-HU  I  )*X2*(CG(  I  1  }  -  I  X1  +  X3I  )*  (Kb 
Jn-K2*vj(  inn _ 

I  ^Gb  +  I  Ft  (  I  I  )  *K  6  %  X  2  *  5*‘  «.  +  H  (  II  )  *  X  2  *  I  K  5+  K  2*  w  (  III  }  ) ^ F  F i / uG  ^ 2 
G  G  —  X  3  +  K  2  t\i  (  I  I  i 

D  X 4  a  i  =  (  k  1 *  n  {  II  I *K6  ^2 * 0X5X 1 )  / GG _ 

DX4  X  2~  K6  *  X  5  /  GG+K  6  *  X  2  *D  X  5  X  2/  G  6 
UXHX;>;=Ko;!!':X2*CxGX3/oG 

UX‘«m  =  1 11  (  II  )*Xl  +  Kt>*X2*bX5Kl)/GG _ _ _ 

DX4K2=  (  K6*X2*LiX5K2*  (  *5  +  K2*.m  I  M  -  w  ti  I  )  *(  K1*H  (  I  l  )  *X  1+K6* 
1  X  2  ^  X  3  I  )  /  o  G  * 

X  *2  _ 

b  X  H  t\  3  =  A  o  *  X  2  *  D  X  3  K  J  /  G  b 
DX^k  +  =  K6*x2*GXGK4/uG 

U  X  H  K  j  -  I  K  o  ;,i  x  2  u  a  z>  k  'j ( is  j  +  K  c  -F  v  »  (  II  )  )  -  t  K 1  *Fi  l  i  i  )  ^  x  1  +  k  o  ^  X  2  5‘c  X  5  ) 
i  I/bb^^L 

LX4ko=X2*I  XJ+Uxd2n6)  /  Go 
0X4K7-K6*X2*GX5K7/GG 
JX { 1 t 1 l < l I )+K5^0X4Xl 
JXI 1  ,2 J=K3^DX4X2 

JX  I  l  >3  J=Kl3» 0X4x3 _ 

JK  I  1  » 1  )  =  -  H  I  I  n*Xl  +  Kb*DX4Kl 
Jk I  I  f 2  J  =  K5*UX4K2 
JK  (  I  ,J>  i  =  Ki*Dx4K3 
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J  M  i  »4)=K3*UX4K4 
oKl  1  »5 )-K5*DX4K5+X4 
JK (  i  = 

J  M  i  » 7  )  =  K35J!UX4K7 

- ^Ai  ^  *•  )  ~K2*U  t  j,l  *i>X  3X  1 +  K  1  *H  {  I  i  )  ^xz^i  1  . 

10+DX4X1+ 

10X3X1) 

JX  (  £,jl2J?K2*W  t  i  i  L*PX4X2+<  K3~K6*X2 )  »£iX5X2~K6:»X5~K7*H  I  I  i  )  * 
I  ft C 0 i 1 l  )- 

11X1+X3+X4+X3)  )  +  K  7  *  H  (  ii  )  £  X  2  #  (  u  A  4X^  +  0X3X*i) 

-  JX  t  2  >  3  )-k  j+wjJj,  )  *uX4A  3  +  (  k  —  K6*  a2)  *0X3X3+k  /*ri  (  1  1  )  »X2*  ( 1  . 
lC  +  u  A^tXi  +  L/ J 
1X3 ) 

A'-t  f,  l  +  {  k  2  -  K  t  *  A  2  )  *DX5K1+  K7*H<  I  Ii  *X2*  (  LX 

1 4K 1  ♦  0  X  5K 1 ) 

JK{  2 , 2  )  =  a  ( l  1  }  *X4+k2*h  (  I  I  )  *uX4k2+  (  K3-K6*X2  )  *k<X3K2+K  7*H  (  I 

_ i  1  )  *>2*(  LX4i\ _ 

12  +  0IX3K2) 

JK(2»3)  =  X3  +  K2*i«{II  )*0X4K3+(k3-K6*X2)*DX5K3+K7*H(  il  j  X  2  ^ 
1 1QX4K3+QX5K  _ 

13) 

II  )*L)X4K4+(K3-K6*x2i5i£LX3K4+K7*H(  I  L  )  *X.2*  l  uX 
_ I4K4+UX3K4J _  _ 

Jklt:  ,  3  I-K^^W  (  II  )^0X4K5+(K3-Ko5i!:X2)*DX3K3  +  K7«Hi  I  1  )  *X2* ( OX 
i4K3  +  JX  3K 3  j 

ji\  {  2  y  t  )  -  A  l *  x  §,  +  t\  a  *  »,  {  i  i  )  *DX4K6+(  K3—K6*X21  *QX  5K6+K  7 *J«| |  j  I )  # 
1X2*  i  oA^Ko  +  U 
1X3KE  ) 

- J.k  (  2  t  7  J  —  H  1  1  I  J*X2*(Cutl  I  )-(Xi+A3+X4+X3)  ) +  K2  *  a (  I  I)  *UX4k  /  + 

1  (  K  3  -  X  o  *  A  2  )  * 

10X3K  /+K7*H(  I  I  j*X2*  (OX4K7 +0X5X7  ) 

JX(3» 1 i«&4*DX5Xl 
ja( 3,2 J=K4*0X5X2 
J  X  (  3  »3)=K4* 0X5X3 

JK  (  3  1 1  )  =  X4» 0  X3K 1 _ __ 

Jk( 3  )=X4*0X3K2 

J  i\  (  3  » 3  )  =  Ah''4uXjK3 
J  X  (  3  »  4  )  =  A  J  +  k  h  *  U  X  3  X 4 
J  K  1  3 , 3  )  -  X  4  *  u  .A  5 X  3 
JK  (  3 r6)=K4*0X5Ko 

Jk  (  3  t  / )  =  K‘-t*U  X3K7 _  _ 

13  RETURN 

END 


1  i  t  1  )  i  L> 

+ 

k  i  ) 

A  —  l  \  t  i  /  .» 

it  I'.  I  )  '  •  '  .  ~  l  .  '  ' 

r  1  «  A  O  *  - 

l  .  .  C  ; 

t  "  ;  '  a  •  >  .  .  I  ■  -  l  f  .  i  . 

~  \  i  »  I  -  I  i 

>  \  '  [  .  ')  1  l  I  l  V  <  »  t  1  ,  . 


V 

V 

i  ■ 

-  t  -  t 

)  k. 

\  „  '  a 

r 

i 

i  / . 

i  > 

.  i 

/I  —  l  t 

) 

i  i,Vv.\g 

■  i 

• 

i 

*.  »  J  •  ■  1  i  | 

\Ai. 

\  ) 

■■  \  .  . 

*  - 

-  ,i  ; 

t 

/ 

-  '  -  V 

/  . 

t  i 

i  > 

r 

\ 

V 

) 

)  w 

t 

f  -.Hi 

t  1  f 

i 

)  * 

V 

t  . 

l 

1 

t  C  /I  A  k.  1 

1 

4 

l 

) 

1  1 

\  i 

.  f 

\ 

L  t  i  i  1 

)  /  o 

'  \  I  .  I  •  \  ,  i  '  i  •••  i  l  /tC  •  i 


X 

V 

X 

i  L 

r 

. 

X 

i  A  o 

\ 

X 

X 

I  ll 

J  •  V- 

- 

t 

f 

- 

)■  \ „ 

t 

/l  * 

~ 

l 

X 

J  •  o 

V 

- 

i  ■  U. 

l 

V 

E  -2  1 


SLdRuuT 1 NE  SLCPt ( K  ,  T  ,  1 oACK ) 

S.LQ.PE  &ji&R.Q,U,,T IHt  EQB.  i  < » ._  >•  QU&  constant  chemoell 
PRO 6k E ,v  AT  80  DEGREES, 

KEAL^d  X(iC,C5),jX{05,05),JMQ5,o5),b(05>, (,C(11,Q5),X0 
_ 1  (  30,05) _ 

ii\iT  EGER  K  *  NE  »NK  ,  NO  ,  1  1 ,  NOS  »  i  JACK 

Cwvi  JN  /  A  t:3  /  C ,  NE  ,  NK  ,  NO /A  35/  Ja  ,  JK,X,G/Ab8/XC,il,NDS 

KE  AL  vd  X  1,A2, X3  , A 4  >X3i  XC  »Ki  ,  K2  jR3*.K4  ,  K3_*.JK6 j.K 7  *  n  1 05 )  ,  x  {  0 

15)  ,  CC { ) 

1 , FF , ub  ,GG2 ,0X4X1 ,0X4X2 ,0X4X3  ,DX4k1  ,DX4K2 ,DX4K3 ,QX4K4,DX 

_ l^K  5 _ 

IN  1  t.GEH  i  ,  j 
DO  16  1  =  1 ,  N  E 

ir  (L.4rtj(  X  (  R  til)  .  L  E  .1  NO)  GOT  l  16 
*RlIEt6»i7* 

STOP 

I  u  C  C  r>.  I  1 0  C  E _ _ 

II  FuRMAT  (  'OsIATE  VARiAoLfc  BQUNu  EXCEEDED’) 

X  1  —  X  (  R  ,  1  ) 

~ - X2Lfr_XJ,ii  1..2.J _ _ _ _  ________ 

X3- X I R  *3 1 
K 1  =  0  ( k  ,  1  ) 

_ K2  =  (.  (  R  ,  2  ) _ _ _ _ 

K3=C(K,3) 

K5  =  G ( R  ,  4 ) 

GOT  C  111,12.)  ,  No  . . . 

11  LO  i 3  1=1, N t 
G( 1 ) =0 .0 

_ DO  14  J  =  1  ,  N  E _ _ _ _ 

14  JX( i , J)=C.0 
DO  13  J  =  1  ,  N  K 

12  jr  (  i ,  j  )  =  c  . ; 

NG=  2 

K4  =  C.C 

_ u*l-\  H  { i  )  ,  a  (  l  )  ,  00  (  1  )  /Q  .  Q  1 2  / ,  5 5 . 6  ,  C  .  0  79  1  / _ 

GAIA  Hi  2 ) ,'rt  ( 2 ) , COl 2)/0 .0230, 56.6 ,0.0795/ 

DaIA  H ( 3  j  ,  to ( 3  )  ,00(3) /0. 0  36 3,  35.6,0.0822/ 

DATA  H (4 ) , rt ( 4 J ,  CO ( 41 1 0 .0345 #55.6,0.0825/ 

DATA  H ( 5  )  ,  W ( 5 ) , CO ( 5 ) / C .0698, 55.6,0.0835/ 
i  2  X5  =  OL(11)-a1-a3 

Ft~  =  Kj*H  (  1 i  )  *X2*X5  +  K1*H U 1 ) *W (  1  I  ) *Xl _ 

G  G=  K  2  +  K  3  +  K  4  *  X  2  +  K  5  *  n  ( 11 )  *X2 
X4=FF / GG 

A  o  =  A  j  -  X  4  ________ 

G  (  i  )  =  “  K  i  *  r  i  l  1  )  ♦  W  (  I  i  )  *X1  +  K4 Al;':a4 

j(e  )  =t\  2  $  X4 +K  1  ^  h  (  11  )*vA  II  )  *X1-K4*X2*X4-n5*H  (  I  1  )  *X2*X6 

0-  (  J  )  =K3»X4 _ _ 

IF {  I  JACK  .01  .0 )  GOTO  15 
002=^3* *2 

0X4X1=  (  K  1  *H_(  1 1  )  l  1  I  )  -  K  5  *  H  (  i  i  )  *X2  )/ub 
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ua^Xc—  l  K  2  $  Fi  (  I  i  i  *X5*GG-  Ik4  +  K3*fi{  i  1  )  )*|-F  )  /ob 2 
&X4X3^K5*H  ill)  ■>' Ac  /  lu 
0X4K1=H( I l ) *W( IX|*X1/GG 
Li  A  4  is  £  =  -  i  •  C  /  u  G  2 

la4nJ=uX^K/ _ 

CX4K4=X2*Ua4K2 

LiX4K5=  ih  (  1  i  )  *X2*X5*GG-H(  I  I  )  *X2*EF)  /  Gb2 

11=^1  tl  i  )  *  W  I  II  )  +K4»X2*iJX4  A  1  _____ 

JX(  1  t  C  J— 'K4  ^  X  <i+KM;A23';UA4A<L 
JX{1  ,3)=K4*X2*DX4X3 

.  JM  1  t  1  J  =  ~h  i  1  1  )  »‘W  I  I  1  )  *X  1  4-K4*X/!*iiX4K  i _ 

JM  i  »2)=K4*X2*Da4K2 
J|\(  1  f  3  )  =K4*X2*0X4K3 
JK-U.  »4)~K4*X2*0X4K4 
FF=?K2~tK4*X2+K5 *  H  C 1 1 )  *X2 

oX( 2  » 1 )  =  K 1 #Fi { II )  I  1 i+K3*H(  Ii  I*X2  +  FF*L)X4Xi 
-J-X  (  C.  ,c  J=- iW*X4-Kb*Ht  i  I  I  frXt^  +  F  F»nX4X2 _ 

JX(2,3)  =  FF*CX4X3  +  Kb*Fi(  i I ) *X2 
J  M  2  >  1  J  =  FF*DX4Kl+rii  1  i  )  *  W  l  L  1)^X1 

J  M  4  ,  l  !  =  r  r  *  u  a  ±  &  2  »X4 _ _ _ _ _ _ _ 

J  Kr  (  2 , 3  j  =  F  F  *  0  X  4K 3 
J  M  2  >4)=FFj!cUX4X3-Fi(  1  I  )  *  X  2  *  X  6 

■  JX(  3  .  1  )=K.^*iJA4a1 _ __ _ 

JX  (  3  , 2 I=k3^DX4a2 
JX{ 3 ,3 ) =K3*DX4X3 

J  M  3  / 1  )  -  ivJ  *U,X4iU _ _ _ _ _ ___________ 

JM  3  »  2  }=*K3*0X4K2 
JK ( 3  » j )■=  K3*DX 4K 3  +  X  4 

JM  3  .4  )  =  K3*uX4i\t> _ _ _ _ 
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S bb K  J 0  I  INC  SLoPf:{k»T,ijAtM 

b  bbu^L  juuiMjUi  iNt  Tpk  Trie  FOUR  CONST  ANT  CH  EMC  ELL 

C  akoplch  AT  100  DEGREES* 

k  t  A  L  *  d  A  i  1 0  »  0  b  )  >  J  X  (  Q  b  » 0  5  )  »  J  K  {  0  b  »  0  b  ]  »  G  (  G  5  )  »  T » C  (  1 1  » 0  b )  »X(J 

_ LIavPjjJ _ _ 

INTEGER  R  »NL  »  NK  »  NO  » 1  I  f  NOS  »  I  JACK 

COMMON  /AEi/b»NE»NK»N0/AB5/JX* JK»X,G/ABb/XC, li ,NOS 

'V  b  A  b  J  A  i  i  A  /_  j  A  j>  f  A  ‘t  j  A  A  »  A  o  J  K  4  »  K  2  J  K  3  J  K  4  »  K  ^  i  K  p  »  K  /  »  M  (  0  b  )  *  W  C  Q 
i  c  J  j  C  0  (  0  6 ) 

i  »  F  P  t  G b  *  b  G  2  »  D  X  a  X  1 »  u X4  a 2  »  0  X 4  X  3  » D  X 4 K  1 » D  X4 K2  »  D  X 4 K  3  t  0 X 4  K4  »  0 X 

_ L  4K  P _ 

INTEGER  ljJ 
pU  Ip  I —  1 »  N  E 

LE  (l)aJp(  bOTb  lo 

WRITE  (*6  *  1 7J 

STOP 

_ lfc  oouTiNbL _ 

1/  H  0  R  R  A  T  ( 1 0  S  T ATE  VAKiAbbE  b  o  U  N  o  LXpElDED*) 

X 1=  X  (  R  »  1 ) 

A  2  —  X  (  K  f  2  J  _ _ _ _ ™™_ 

A  P  -  A  {  i<  ,  J  I 

K 1  —  b  (  R  » 1  J 

_ K  2  =  b  {  R  »  2  ) _ 

Kb  =  C ( K  »  3  i 
X5=b (k,4) 

_ GOTO  ( 11?»  12)  » NO _ 

11  p  0  1 3  I— 1 »N£ 

G ( 1 ) =0  .C 

_ CG  14  J  =  1 »  N  E _ 

14  JX (  I  f  J  J  =0 . 0 
00  13  J  =  1 1 N  K 
13  JM  1 ,  J  )=C*Q 
NO- 2 
K4=  G • 0 

_ GAIA  H  U  i  >  w  tl)  y COil) /Q  .0  12  y  bb.  o  >  0  .  C  7  39/ _ 

DATA  Hi  2  )  ( 2) »C0( 2)20.00241 , 5 5.  by  0.07  74/ 

DATA  H  {  3  )  » w l 3 )  »bGl3J/0*QQbl4»bb.6yC.Q828/ 

DAT  A  h i 4  }  »  a ( 4 ) »CO( 4i /0  *00986  >  55  <  6  y J . 0  7  8  4 / 

PA  TA  H< 5 )  « W ( 51 ,CG(5)/0.0125,55.6, 0*0 749/ 

UaTA  H  i  6  )  i  ,v  (  o  )  »CUl6J/0.0103i55.6fC*013t)3/ 

_ lb  A3  —  Col  i  1  )  ~ A  1  - X p _ _ 

FF=Rb i 1 i )*x2*Xb-*-Kl*Ht 1  I ) ( 1 1 ) *X1 
Gu  =  K2+k3  +  K4,!<  X2  +  Kb*H  (  il  )  &  X  2 
a  4  =  F  F  /  GG _ _ _ _____ 


Ab  =  A  3- >4 

ui  1  )  -  -  k  1  *  h  l i 1 J  *  to (  I  I  ) 
G  {  2  1  -  is  2  *  X  4  t-  m  *  H  i  11  )  * 

-•‘A  i  +  K4*X25I;X4 

to  i  ll  j  *  X  1  -  K  H  *  X  A  * X 4 -  K  b  *H  (  I  i  )  *  X  2  #  X  u 

b  (  3 ) -K  3  *  X4 

It-  t  I  JAOK  .  L  i  .0  )  GOTO 

lb 

t  '  *  )  X  ' 


.  i  ■  t  I  '  U  I  V  •  -  /  ''  -i  1  '  ' 

»  c  /  i 

•  \  \  I  w  \  V  t  t  *  ■>  \  \ 

i  I  \  i  t  V.  I  ■  I  f  t 

\  I  J  t  l  i 

t  .  -  t  i  >'■ 

t  i  -  i 

l  •  .  .  ,  *  l  I  X  i  I  :  )  •  I 

l  I  1  (  I  1  I 

, 

l  I  i 


t  *  C  C  f  -  i  >  v  (  I 

\  \  \  .  1  .  •.  ?  1  .  \  l  i  ^  \ 

\  ,  f  .  ^  !  '  .  \  i  I  l 

\  •  .  \  [  r  1 

'  1  ;  i  .  .  \  \  ‘  J  « 

i  •  I  t 


•  K  K  i  ) 


\  I  '  , 


» 

) 

1 

t 

1  .. 

( 

,  X 

- 

t 

1 

t 

t 

1 

i  . 

V 

f 

i  ~ 

t 

t  l 

f 

• 

-  1 

i  /  .» 

_ 

• 

U  t  i 

< 

• 

0  t  > 

i  .W 

.  .  j 

) 

t 

) 

) 

l 

l 

) 

t 

V 

1 

i  < 

1 

t 

t 

L 

>  -  . 

l 

\ 

- 

--  l 

.  1  J 

V 

1  ■ 

-  V 

• 

I  w1  l  l 


t  -24 


L,X4Al=  (Kl*h(  1 1  )  (  I  I  )-K5*H(  i  I  J  *X2  J/GG 

i  K34ril  i  i  I  &X54GG-  tK4+K5*H<  i  i  /  )*FF  I/GG2 
GX4X3=~K5*H  (  i  i  J  -  A  2  /  bo 
DX4Ki  =  M  il)*k(ll)*Xl/bo 
UX4K^=-i  .Q/ub2 


C  A  4  A  3  = 
L  X  4  K  4  = 
DX4K£= 
JA  (  1 
JA  (  1 

JXt  1 


JA  i  ij  1  >  =  —  H  (  I  I  )  *W  (  I  i  )  *X  1+K44X24UX4K  1 
J  A  1  i  >2i  =  A4*A2*L;X4K2 
jin  i  i  ,  j  )=K4*X2*DX4K3 
JK(  1  »  4  )  —  a  4 a  2 U  X  4  a  3 
H-=K4-a44X2+K54H (  I  1)4x2 

J  A  t  2  >  1  )  —  A  i  4  H  (  1  1  )  4  yy  (  1  I  )  +  i\  3  4  H  (  11  )  *  a  A  +  f  F4jJ  X  4  X  I 


JA  ( 
JX  ( 

JK  ( 

JA  ( 
JA  l 
2K L 


4 

2 

2 

2 

2 

2 


u  A4K2 
X  2  4  L  X  4 
(Hj  i  i  j 
)  —  —  K  1 4 
)  =  A  4  4  X 
)  =  K  4  »  X 


K  2 
4  X_ 

h  t 

-t  + 

2  4 


2  v  a  ;>  o  o - 1 1 1  I  i  )  a  2  4  £  F )  / GG2 
i  i  )  4  > ,  {  i  I  )  +  K  4  A  2  4  0 X  4  A  1 
K4*X^4bX4X2 
0  X  4  X  j 


2)  =-K44 A4-K54H( I i J *Xo+PF*UX4X2 

3)  =  Fi~4QX4X3  +  K34Hi  li  )  4X2 

li = e t * \x a 4a  i + r  (  n  >  4 Vv  { i  n*xi 

2 } =Fr 4DX4K2+X4 

3  )~Ff-4iJX4A3 

4  J  =  Fr»U  A4K3-H  (  1  fi  *X?»Xtj 


JXt  3 

i  )  =  K  3  4  D  X  4  a  1 

JA  (  3 

2  )  —  K  J  4  u  A  4  X  2 

JX(  3 

J  )  =  K  3  4y  a  4  X  3 

JA  t  3 

1  )  =  A  3  4  U  a  4  K  1 

JK  i  3 

Zj=K34uX4K2 

J  K  {  3 

3  )  —  K  J  4  Li  X  4  K  3 

JK  {  3 

4  )  =K34C X4K5 

13  He  1  L  K  N 

END 
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APPENDIX  F 

Tracer  Problem 


In  order  to  correlated  the  behavior  of  the  flow 
of  a  tracer  through  a  packed  bed,  the  following  mathe¬ 
matical  model  was  used. 


dx^ 

dt 


x 


2 


dx2 

dt 


alX2  -  a2  LX1  "  f<t  ‘  V] 


(F-l) 


The  state  variable  x^  is  the  outlet  concentration  of  the 
tracer.  The  inlet  concentration  at  a  time  t  is  given  by 
f (t) .  The  time  delay  tD  is  not  known  but  must  be  assumed 
in  order  to  use  the  proposed  method.  The  optimum  time 
delay  was  determined  by  search. 

To  solve  for  the  parameters  a^  and  a^  was  not  diffi¬ 
cult  once  a  time  delay  was  specified. 

Setting  the  maximum  step  size  S  to  0 . 1  and  the 
^  c  max 

initial  guess  to 


a 


1.0 

1.0 


(F-2) 


■ 


F-2 


the  solution  converged  rapidly,  as  shown  in  Table  F-l, 
which  gives  the  results  for  the  optimum  time  delay 
(tD  =  0.187)  . 


TABLE  F-l 


Tracer  Problem  tD  =  0.187 


Phase  I 
Step  No . 


Sum  of 
a0  Errors 

Squared 


0  1.0  1.0  0.116 

1  5.47  5.24  0.032 


Phase  II 
Iteration 


0 

5.47 

5.24 

1 

9.42 

11.43 

2 

6.34 

10.31 

3 

6 . 54 

10.62 

0.00035 

For  the  other  values  of  the  time  delay  consider¬ 
ed,  rapid  convergence  was  also  obtained. 

The  SLOPE  subroutine  required  for  this  problem 
is  listed  on  the  following  pages. 


' 


F  ~  3 


SUbKGUT  Ii\i£  l  L  >_P  1 1  R  ,  f  ,  I jAlk ) 

SLlfL  jUdklU  r  I  iMc  Fur*  7  He  TRACER  F  RUBLE*!. 
kEAL  XilC»10)»JX(  10,lG)»JK{lC,10ltGilQ),T,L(  11,10) , X 0 ( 3 
1C, 1C) 

- INTEGER  R,N£ «NK.NQ. I ; T  IDS . I JACK _ 

COMMON  /AB3/C  ,  NE , NK ,NG/ AB5/ JX , JK » X ,G/Abb/XC , 1 1 ,  NDS 
REAL  X1,X2,K1,K2,K3,K4,A3,X4 
I N 1  E  u  E  R  i,J 

u  L  lc  i=l, kE 

IMABS(X(Rr  I  )  )  .Lfc.  10.0)  GGTu  16 

- IolEJCLl:  (6,17) _ _ _ _ _ 

STOP 

10  CUM  IfsiUE 

1  t  FORMAT  {  'CGTATc  VAkiAoLt  BOUND  EXCEEutU*) 

X1=X ( R  ,  1 ) 
x^=xik,2) 

kl=C(k,l) _ 

K2  =  C l R  »  2 ) 

GO  TU ( 1 1 , 12 )  ,  NO 

11  Cl  i 3  1  -  1 ,  n  E 
G(IJ=0.0 

Cl  1a  j  —  1 ,  N  E 

1 A  jX(l,J)=C.j _ 

Cu  13  J  =  1 ,  N  K 
13  JK i 1 , J ) =0 . 0 
Nu  =  c 
A  =  C . 186 

12  311  =  31 CT-A  ,  I  I) 

_ CJ  1  )  ~  X  2 _ _ 

l ( 2 ) =-Kl *X2-K2* ( X  1-3 1 i ) 

1 F {  I  JAlK  .  l T  .  C )  RETURN 
J  X  (  1 ,  A  )  =  1 . 0 
J  X  (  c  ,  1  ) =-R  2 
jX { 2 , 2 ) =-k 1 

_ JK(2tl)  =  -A2 _  _ 

JK  ( <£  ,2)=-Xi  +  Sll 

RETURN 

END 
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REAL  huNCTIGK  S1(T,1) 

L  LmPLI  t-  UJSi  l  T  I  Q.&  REwUikEb  bV  The  SLUPt  SUBROUTINE  PGR 

c  the  tracer  problem 

S 1  =  C .  7  5 

_ IF(T.lt.q.o)  ij o t n  in _ 

GOTO  111 ,  12, 13, 14) , I 
1  i  al=Sl+C.25*slN(10.0*T) 

GOTO  IC 

12  Gl  =  El  +  C.23^ElNl 1 . C  *  T  J 
GUTG  1C 

- U  Sl  =  SlfrC.25*SlNtG.l  *T) _ 

GOTO  1C 

14  S1=S1+C.25*SIN( 5.0*T) 

LQ  RcTCRN 
END 


I  ,  i  i. 

i  •  •  .  u  .i 
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.  1  /  1  C  .  +  i  <  I  <i  i  - 

l  ♦  i  /  1 1.  ■  .  •  .  - ,  u 

.  . 

I 


